INERALIZED MOMENT METHODS 
ELECTROMAGNETICS 

rmulation and Computer Solution of Integral Equations 
anson J. H. Wang 


is comprehensive treatise covers the broad class of computational 
‘hniques known as the method of moments, techniques widely 
ed in calculations of electromagnetic problems. The book presents 
t only a unified approach to the numerical aspects of moment 
+thods, but also the analytical background, including an up-to-date 
dcomprehensive look at electromagnetic theory. Designed as both 
rofessional reference and graduate level textbook, the book’s prac- 
al, hands-on format includes exercises, four computer programs 
mplete with user’s guides, and examples with solutions. 

91 (0 471-51443-8) 553 pp. 


EECTROMAGNETIC WAVE THEORY 
cond Edition 
1 Au Kong 


.. presents a unified macroscopic theory of electromagnetic waves 
topics essential to the understanding of electromagnetic waves are 
io selected and presented. ..a useful graduate text.” 

—Physics Briefs 
‘cognized for its modern, unified approach to the subject, this text 
s been used by countless first year graduate students as an intro- 
ictory reference, With an emphasis on mathematics, problem solv- 
g, and physical interpretation, the book’s wealth of examples 
vers guidance, propagation, and radiation of electromagnetic 
aves; metallic and dielectric wave guides; resonators, antennas, 
id radiating structures; and more. Organized along lines of increas- 
g technical and conceptual complexity, the book’s format allows for 
dependent teaching of topics. 
90° (0 471-52214-7) 704 pp. 


UMERICAL TECHNIQUES FOR MICROWAVE 
ND MILLIMETER-WAVE PASSIVE STRUCTURES 
lited by Tatsuo Itoh 


(his book will go a long way to help introduce the engineers of this 
eration (and the next generation) to a valuable set of tools.”’ 
—IEEE Antennas and Propagation Society Newsletter 


\is unique compendium of the most widely used numerical meth- 
ls for analyzing the passive structures appearing in microwave and 
illimeter-wave integrated circuits represents the collective expertise 
an international team of engineers. The book includes useful tips 
1 writing computer programs based on the method described and 
cludes sample computer programs for reference. 

89 (0 471-62563-9) 707 pp. 


ISBN O-471-528D4-8 


WT 


ILEY-INTERSCIENCE 

hn Wiley & Sons, Inc. 

ofessional, Reference and Trade Grou; 

Third Avenue, New York, N.Y. 10158-0012 


Maveb «Chichacter « Rrichane s Taranta 


9014 


nN 


JOUSPLUOJ) 


— 
=. = 
ee 
f=) 


SOABMOJONA 2!” 


jndwio 


2 


‘SPOUVN | 


Computational 
Methods for 
Electromagnetics 
and Microwaves 


Richard C. Booton, Jr. 


) OPTICAL ENGINEERI 


the WILEY SERIES IN MICROWAVE 
KAI( as 


mputational Methods for 
omagnetics and Microwaves 


= 


en i 


WILEY SERIES IN MICROWAVE AND OPTICAL ENGINEERING Co m p utati O n al 

rind Methods for 
Electromagnetics 

maar TO ELECTROMAGNETIC COMPATIBILITY and M i C rowaves 


OPTICAL COMPUTING; AN INTRODUCTION 
Mohammad A. Karim and Abdul Abad S. Awwal 
COMPUTATIONAL METHODS FOR ELECTROMAGNETICS 
AND MICROWAVES 
Richard C. Booton, Jr. 
FIBER-OPTIC COMMUNICATION SYSTEMS 
Govind P. Agrawal RICHARD C. BOOTON, JR. 
Department of Electrical and Computer Engineering 
The University of Colorado at Boulder 
Boulder, Colorado 


®) 


A WILEY-INTERSCIENCE PUBLICATION 
JOHN WILEY & SONS, INC. 
NEW YORK / CHICHESTER / BRISBANE / TORONTO / SINGAPORE 


> 


In recognition of the importance of preserving what has been 
written, it is a policy of John Wiley & Sons, Inc., to have books 
of enduring value published in the United States printed on 
acid-free paper, and we exert our best efforts to that end. 


Copyright © 1992 by John Wiley & Sons, Inc. 
All rights reserved. Published simultaneously in Canada, 


Reproduction or translation of any part of this work 
beyond that permitted by Section 107 or 108 of the 

1976 United States Copyright Act without the permission 
of the copyright owner is unlawful. Requests for 
permission or further information should be addressed to 
the Permissions Department, John Wiley & Sons, Inc. 


Library of Congress Cataloging in Publication Data: 
Booton, Richard C., 1926- 

Computational methods for electromagnetics and microwaves/Richard C. 
Booton, Jr. 

Pp. em. —(Wiley series in microwave and optical engineering) 

“A Wiley-Interscience publication.” 

Includes bibliographical references. 

ISBN ()-471-52804-8 

1, Electromagnetism—Data Processing. 2. Microwaves—Data 
processing. I. Title. Il. Series. 


QC760.54B66 1992 91-34450 
621.3—de20 cIP 


Printed in the United States of America 
198765432 


1.1. Electromagnetic Problems Considered 
in This Book 

1.2. Basic Numerical Methods 

1.3. Solution of Algebraic Equations 

1.4. Accuracy Considerations and Richardson 
Extrapolation 

1.5. Integration Example 

r Problems 

Computer Project 1-1 

References 


Difference Method 

2.1. Finite Differences in One Dimension 
2.2. A One-Dimensional Differential Equation 
__ Example 

2.3. Finite Differences in Two Dimensions 
24. Two-Dimensional Capacitance Example 
2.5. Open Regions 

2.6. Generalizations 

Problems 


vi CONTENTS CONTENTS 


Computer Project 2-1 38 7. Nodal Expansions and the Weak Formulation 
References 39 Time-Dependent Variables 
Problems 
3. Finite-Difference Determination of Eigenvalues 41 Computer Project 6-1 
. : . ‘ References 
3.1. Eigenvalues in One Dimension 41 
3.2. Waveguide-Mode Example 42 Mom 
3.3. Numerical Evaluation of the Determinant 46 ad sai 
3.4. Iterative Solution Methods 52 7.1. Linear Operators a. , 
Problems 56 Approximation by Expansion in Basis Functions 
References 58 Determination of the Parameters 
Differential Operators 
sie ri " Integral Operators 
4. Finite-Difference Time-Domain Method 59 Pulse Functions 
4.1. Wave Equation in One Spatial Dimension 59 7. Parallel-Plate Capacitor in Two Dimensions 
4.2. Time Quantization 60 Problems 
4.3. Initial Conditions 63 Computer Project 7-1 
4.4, Waves in Two and Three Spatial Dimensions 64 Computer Project 7-2 
4.5, Maxwell’s Equations 68 Reference 
Problem 72 
roe in cs aaa Fo ering Solutions by Method of Moments 
Fundamental Scattering Equations 
5. Variational and Related Methods 75 Scattering of a Plane Wave off a Thin Cylinder 
Scattering of a Plane Wave off a Thick Cylinder 
5.1. Stationary Conditions for Functionals 75 Problems 
5.2. Example of Electrostatic Field Energy 78 Computer Project 8-1 
5.3. Variational Expression for Eigenvalues 79 References 
5.4. Rayleigh-Ritz Method 81 
Problems 83 tral Analysis with Fourier Series and Fourier 
Reference 85 q 
525 A. Basic Fourier Series Relations 
6. Finite-Element Method 87 Example Involving Laplace’s Equation 
6.1. Basic Concept of Finite Elements 87 Inductive Iris in a Waveguide 
6.2. Finite Elements in One Dimension 89 1.4. Basic Fourier Integral Relations 
6.3, Linear Interpolation for Isoceles Right Triangles 92 .5. Fourier Integral Solution 
6.4. Square Elements 98 Problem 
6.5. General Triangular Elements 100 Computer Project 9-1 


6.6, Higher-Order Interpolation with Triangles 104 References 


viii 


12. 


10. 


11. 


CONTENTS 


Spectral Analysis of Microstrip Transmission Lines 


10.1. 
10.2. 


10.3, 


10.4. 


Microstrip Transmission Line 

Quasi-TEM Spectral Analysis of Shielded 
Microstrip Transmission Line 

Quasi-TEM Spectral Analysis of Unshielded 
Microstrip Transmission Line 

Comparison of Fourier Series and Fourier Integral 
Solutions 

Problems 

References 


Spectral Analysis of Microstrip Circuits 


11.1. 
11.2. 


ILS 


11.4, 


Microstrip Circuits 

Quasi-TEM Spectral Analysis of Shielded 
Microstrip Circuits 

Quasi-TEM Spectral Analysis of Unshielded 
Microstrip Circuits 

Full-Wave Solution of Unshielded Microwave 
Circuits 


Mode Matching 


12.1. 
12,2, 
12,3; 


Approximating a Function by a Series 

Inductive Iris in a Waveguide (Revisited) 

Other Examples of a Waveguide Mode Matching 
Problems 

Computer Project 12-1 


Concluding Comments 


Index 


179 


181 


Preface 


This book is based on my experience teaching a beginning graduate 
course in the Electrical Engineering Department at the University 
of California, Los Angeles, and in the Electrical and Computer 
Engineering Department at the University of Colorado at Boulder 
and is an expansion of notes prepared for this course. This book is 
intended to serve as a textbook, not as a reference. My experience 
is that there are several excellent books that can be used as 
references and assigned as supplemental reading, but that a real 
need exists for a textbook that can be used for an electrical 
engineering course. I have deliberately kept the length and depth of 
the coverage to what can be covered in a one-semester course. 

I firmly believe that numerical methods can be learned effec- 
tively only by solving realistic problems, and such solutions require 
that the student write and run computer programs. Although useful 
Programs can be written in BASIC, PASCAL, and other languages, 
almost all electrical engineers use FORTRAN, and examples in this 
book use Microsoft FORTRAN. Although the question of what 
computer language ought to be taught and used is a controversial 
Subject with many people, it seems to me not to matter much what 
language is used. Far more important is the understanding by the 
Student of the algorithms and how they relate to the physical 
problems. In a single university quarter or semester, three or four 
good computer projects can be carried out. I have included what I 
think are relevant project descriptions in the appropriate chapters. 
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The enterprising instructor can easily create others. The essential 
point is that the student must write and use a workable program. 
Although later in his career he probably will use professionally 
written “canned” programs, the emphasis here is on his under- 
standing the basic algorithms. 

The emphasis here is on electromagnetic and microwave prob- 
lems and the fundamental algorithms that can be used as the basis 
for computer programs that produce useful numerical results. Many 
such programs can be run on personal computers. Of course, more 
memory and speed are needed for complex problems. I spend little 
time with students on software niceties, and the better students 
seem to pick up such capabilities elsewhere. 

I wish to express my deep appreciation to four people: Zdenék 
Kopal, who first taught me the fundamentals of numerical methods 
many years ago at MIT; David Chang, who reawakened my interest 
in numerical methods for electromagnetics; Nicolaos Alexopoulos, 
who encouraged me to create and teach the course at UCLA; and 
especially to my wife, Patricia, both for her general encouragement 
and for her toleration of my many hours at the computer in our 
bedroom. Also, I wish to thank Kai Chang for his decision to 
include this book in his series and for his helpful suggestions on 
material to be included. 


Ricnarp C. Booron, Jr. 


CHAPTER ONE 
—SSa 


Introduction to Numerical 
Methods 


1.1. ELECTROMAGNETIC PROBLEMS CONSIDERED 
IN THIS BOOK 


The numerical solution of partial differential equations and of 
integral equations has a long history and a variety of methods have 
been developed. Electromagnetics leads to such equations, and 
numerical studies in electromagnetics have been responsible for 
much of the current development of improved and more powerful 
methods. The three fundamental methods of finite differences, 
finite elements, and moment methods are concentrated on here, 
with brief coverage of some other techniques. Emphasis is on the 
electromagnetic problems and the algorithms that can be used as a 
basis for computer programs that yield useful numerical results. 
The programs described in this book are written in Microsoft 
FORTRAN and have been run on several personal computers 
using the Microsoft MS-DOS operating system. Other languages, 
such as BASIC and PASCAL, can be used, although special modi- 
fications will be required when complex variables are involved. With 
appropriate modifications if different FORTRAN versions are used, 
the programs can be run on almost any personal computer, work- 
Station, or other computer. Adequate memory and processor speed 
are required to achieve high accuracy. Real physical problems 
involve three spatial variables and one time variable. Because the 
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full treatment of such problems uses large memory and long run 
times, wherever possible simplification to one or two spatial vari- 
ables is carried out. 

The principal engineering applications of numerical methods to 
electromagnetics are to guided waves, antennas, and scattering. 
Guided waves probably offer the easiest introduction to most of the 
numerical methods. The analysis of TEM waves on two-conductor 
transmission lines leads to static calculations of capacitance, which 
are used to illustrate the three basic methods. Analysis of wave- 
guides leads to eigenvalue problems which we will solve with both a 
classical method and a newer method. The first uses algebraic 
eigenvalue techniques coupled with the finite-difference method, 
and the second approach uses the time-domain finite-difference 
method coupled with the discrete Fourier transform. Except for the 
latter method, these methods are classical methods with long histo- 
ries. 

The analysis of microstrip and similar transmission lines, which 
are becoming increasingly important in microwave circuits, is a 
more difficult problem. Because microstrip and similar lines involve 
nonuniform dielectrics, they cannot support a TEM wave. At low 
frequencies, the dominant mode can be approximated by a TEM 
wave. This approximation leads to a method of analysis referred to 
as quasi-TEM or quasi-static, which is essentially the same as the 
static analysis of two-conductor lines with a uniform dielectric, At 
higher frequencies, this approximation breaks down and more so- 
phisticated methods are needed. Accurate solution of such prob- 
lems results from the combination of the moment method with 
spectral (Fourier series and Fourier transform) techniques. In the 
last few years, a number of other methods have been used for such 
problems. 

Another major class of electromagnetic problems involves analy- 
sis of scattering of waves. The simplest scattering problem involves 
a plane wave in free space striking a body, frequently taken as a 
perfect conductor. The incident wave causes currents in the body 
and these currents create a second wave called the scattered wave. 
We shall illustrate scattering calculations for a problem of this type. 
The analysis of the effects of an obstacle in a waveguide offers 
another example of scattering. 
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4.2. BASIC NUMERICAL METHODS 


“The reader is assumed to have a basic background in numerical 
" methods, such as is provided in the undergraduate course at many 
“universities and covered by a number of excellent texts [1]-[3]. The 
~ most important topic in this background is the solution of systems 
of linear algebraic equations. More advanced topics can be found in 
“several numerical texts [4]-[5]. ; 
Electromagnetic problems are sometimes described by differen- 
tial equations, sometimes by integral integrations, and sometimes by 
‘minimization of an integral such as the energy integral. The un- 
known function is usually continuous and depends on continuous 
independent variables. Computers clearly can handle only a finite 
‘set of numbers, whereas the differential and integral equations that 
~ describe electromagnetic fields, or indeed almost any physical vari- 
able, involve an infinite set of numbers. By one means or another, 
the more accurate continuous equations must be manipulated to 
_ produce equations that involve a finite set of numbers. The final 
equations usually are algebraic and in electromagnetics, usually 
linear, and thus solution of a set of linear algebraic equations is 
involved. Although in this book we consider only situations with 
passive components, which lead to linear algebraic equations, if 
active devices such as transistors are involved, nonlinear equations 
result. Solution methods must then be modified, almost always 
_ through use of iterative methods. 
Use of any of the numerical methods considered here involves 
three steps: 


i 
1 


1. Preprocessing, to derive the coefficients in the algebraic equa- 
tions 

2. Solution of the algebraic equations 

3. Interpretation of the results 


The various numerical methods involve different amounts of theo- 
Tetical effort and different computational efficiencies. As we shall 
_ See, increasing the efficiency of the computational process usually 
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FIGURE 1-1. Two-dimensional potential problem. 


must be paid for with more preprocessing effort and more use of 
electromagnetic theory, 

Although, as indicated earlier, a variety of such methods have 
been developed, three seem to be the most fundamental: the 
finite-difference method, the finite-element method, and the mo- 
ment method. As an example, consider the two-dimensional poten- 
tial problem illustrated in Fig. 1-1. The potential d(x, y) that 
satisfies Laplace’s equation and certain boundary conditions on the 
two conductors is to be determined numerically. The finite-dif- 
ference method concentrates on a finite mesh of points arranged in 
a rectangular grid or mesh, as shown in Fig. 1-2. By techniques to 
be discussed in Chapter 2, Laplace’s equation is used to establish a 
set of linear algebraic equations whose solution gives approximate 
values for the potential on the mesh points. The closely related 
finite-element method divides the region of interest into a number of 
subregions (the finite elements that give the method its name), as 
illustrated in Fig. 1-3. The field in each element is approximated by 
a simple algebraic expression, usually a low-degree polynomial, and 
by techniques to be discussed later, the values of the field on a 


FIGURE 1-2. Finite-difference mesh. 
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FIGURE 1-3. Triangular finite elements. 


finite set of nodal points are found as the solution of a set of linear 
algebraic equations. 

The third method is the moment method, sometimes referred to 
as the method of moments. This method involves more sophisti- 
cated mathematics and in return achieves greater computational 
efficiency. For the potential problem illustrated here, the potential 
is expressed as an integral over the center conductor of an inte- 
grand, which is the product of the charge density and a function 
known as the Green's function. Much of the effort in using the 
moment method involves calculation of the Green’s function, which 
frequently is determined as the sum of an infinite series, The charge 
density is expressed as a linear combination of a set of known 
functions (known as basis or expansion functions). By methods 
discussed in the following section, a set of linear algebraic equa- 
tions are determined whose solution is the set of coefficients in the 
expansion. 

Many other methods, some of which are variations on these basic 
methods, have been developed. Although we devote one chapter to 
consideration of one additional method, even a much larger book 
would be hard-pressed to cover all the methods that have been 
utilized. More detail on the methods covered here and on other 
methods can be found in several books [6]-[11]. 

With any numerical method, maximum use should always be 
taken of symmetry to reduce memory requirements and computing 
time. For example, the two-dimensional potential problem illus- 
trated in Fig. 1-1 has one axis of symmetry. Therefore, for finite- 
difference and finite-element methods, one need consider only 
one-half of the region, as illustrated in Fig. 1-4. At the artificial 
boundary caused by the cut along the axis of symmetry, appropriate 
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FIGURE 1-4, Use of symmetry for potential problem. 


FIGURE 1-5. Problem with two degrees of symmetry. 


boundary conditions must be imposed. Moment and other methods 
also can take advantage of the symmetry, by means to be discussed 
later. Higher degrees of symmetry may exist. For example, Fig. 1-5 
illustrates a situation with two axes of symmetry. If the region is cut 
along each axis of symmetry, only one-fourth of the region need be 
considered, Although rare, sometimes even more symmetry exists. 
In Chapter 2, a situation involving three axes of symmetry will be 
analyzed. 


1.3, SOLUTION OF ALGEBRAIC EQUATIONS 


The systems of linear equations to be encountered in the solution 
of electromagnetic problems fall into two classes. In the first, most 
of the possible terms are zero and the associated matrix is de- 
scribed as sparse. Such systems frequently are solved by iterative 
methods such as Gauss-Seidel [12], but powerful methods for 
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direct solution exist [13]. As we shall see when considering discretiz- 
ing methods such as finite differences, frequently the matrix terms 
are not initially computed and stored but are computed when 
required. This is an advantage when dealing with very large matri- 
a methods, in particular the moment method, involve matri- 
ces that in genera] have no nonzero terms and are referred to as 
dense matrices. Although indirect methods are sometimes used with 
dense matrices, direct methods are more commonly used. Many 
algorithms exist for the direct solution of such equations, many of 
which are variations of Gaussian elimination. All useful algorithms 
utilize pivoting, which involves interchange of rows or columns to 
avoid division by zero or by a small number. Partial pivoting as 
usually implemented utilizes interchange of rows, and full pivoting 
interchanges rows and columns. Partial pivoting algorithms are 
usually adequate and are commonly used, Sometimes the matrix is 
nearly singular, referred to as ill-conditioned, and full-pivoting 
algorithms are preferred. Most computing networks permit access 
to excellent matrix inversion routines, and for stand-alone comput- 
ers without such routines available, excellent routines have been 
described [14]. 


1.4. ACCURACY CONSIDERATIONS AND RICHARDSON 
EXTRAPOLATION 


As discussed earlier, we typically convert the original differential or 
integral equation into a finite set of algebraic equations either by 
concentrating on a finite set of nodal values or by expanding the 
unknown function in a finite sum. At least until we have to worry 
about round-off error, we expect the results we calculate to become 
more accurate as the number N of unknowns increases. Conver- 
gence may be slow, however, and the difference between one value 
of N and the next is usually small. An approach due to Richardson 
frequently is useful. This extrapolation method is described in many 
numerical texts in connection with what is called Romberg integra- 
tion but is a very broadly applicable technique [15]. 

4 One approach to using Richardson extrapolation is to succes- 
Sively halve the separation between points or double the number of 
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functions. It is reasonable to assume that results calculated with 
one level of quantization are at least as accurate as the extent to 
which the results calculated agree with the results at the previous 
level of quantization. We would like to continue in this manner, 
successively doubling the number of points or functions until we 
achieve the accuracy we desire. However, memory requirements 
and computing time increase rapidly as the number of points or 
functions increase, and we may not be able to achieve the accuracy 
we desire in this manner except for simple problems. Consider the 
Ng of a series of finite-difference calculations in the following 
table. 


h Calculation 
1 70,83350 
05 59.87117 
0.25 55.67812 
0.125 53.79773 


Clearly, we could continue to decrease A because the results of 
the calculation have not converged to an accurate enough result. 
Notice, however that there is a clear trend in the results, which is 
made even clearer in Fig. 1-6. Extrapolation to h = 0 is suggested 


0 0.25 05 0.75 10 


FIGURE 1-6. Graph of calculated results. 
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by the figure. A numerical extrapolation developed by Richardson 
+t upon a polynomial fit to the calculated numbers will now be 
explained. Richardson extrapolation has a wide range of application 
a numerical analysis. 

‘As a first approximation, it is reasonable to assume that the error 
for the smallest values of / is proportional to h, which means that 
the calculated values satisfy 


C(h) =a + bh. (1-1) 


Tf calculations are carried out for the two smallest values of h, s, 
and 2s, then 


at+bs=C(s) (12) 
a + 2bs = C(2s), 
from which, by simple algebra, we can compute 
a = 2C(s) — C(2s). (1-3) 


The value of a is, of course, the value of the polynomial for h equal 
to zero, and hence the desired result. This extrapolation uses the 
last two computed values. Since we have more than two computed 
values, we should be able to do better by using a higher-degree 
polynomial that agrees with the computed values. In the calculation 
above we used two computed values and a first-degree polynomial. 
In general, the degree of the polynomial that we can use is one less 
than the number of computed values that we utilize. 

For three and more points, formulas similar to (1-3) can be 
derived, One fits a polynomial to the computed results, calculates 
the coefficients in that polynomial, and then evaluates the polyno- 
‘mial at zero. A more systematic approach can be used, however. 
Suppose that the calculated result is expressed as a power series in 
A in the form 


C(A) = Eo a,h", (1-4) 
a0 
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A first extrapolation can be defined by the formula 
E\(h) = 2C(h) — C(2h). (1-5) 


Substitution of this expression into (1-4) gives 


E\(h) = ey a,[2h" ~ (2h)"], (1-6) 


n=0 


which, because the coefficient of the first power of h is zero, can be 
written 


E(h)=a,+ ¥ a,[2h" — (2h)"]. (1-7) 


n=2 


Note that the leading term in the error is a multiple of A. If we 
apply a similar process to extrapolation of E 1 the result is a second 
extrapolation defined by 


4E\(h) — E,(2h) 


E,(h) = 


(1-8) 


If we Substitute the series into this expression, we will find that the 
coefficient of the second power of f is now zero. We can continue 
this process and define the general extrapolation as 


2"E,-(h) — E,_\(2h) 


E,(h) = —— 


(1-9) 


The number of levels of extrapolation is clearly limited by the 
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jumber of calculated values we have. Applications of this process 
the values in the preceding table leads to the results shown in the 
wing table. 


Extrapolations 
h Calculation E, E; Ey 
1 70.83350 
05 59.87117  48.90884 
0.25 55.67812  51.48507 -$2.34381 
0.125 53.79773 ~—51.91734 —52.06143 —52.02109 


ne would probably conclude from these results that the true value 
the solution is approximately 52.0 with a probable error of 0.1. 
ore extensive calculations on this problem indicate that the solu- 
mn to three significant figures is indeed 52.0. 

This extrapolation process can be inserted into the computer 
program that computes the original results. One forms a matrix, say 
O(i, j), and inserts the calculated results in the first column. The 
r extrapolated results can then be calculated by a simple 
rithm, which in FORTRAN can be written as follows: 


QCi,j=(C*QCi,j-1)-QC1-1,j-1))/ 00-1) 
end do 
end do 
] With any set of calculated results the most accurate extrapolated 
‘alues are the diagonal values Q(j, j) and for the example above, 
value to use is O(4,4), The approach described above in effect 
S a polynomial through the calculated results and evaluates 
this polynomial at h = 0. The same result could be obtained by 
-€xplicitly deriving the polynomial through use of Lagrange interpo- 
lation formulas. 
\ In some special applications, such as integration with the trape- 
-2oidal method the error series can be shown to contain only even 
_ Powers of # and the algorithm should be suitably modified to take 
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advantage of this knowledge. The resulting method is known as 
Romberg integration. Usually, however, we have no such knowledge 
about the error and all powers of fh should be included. 

1.5. INTEGRATION EXAMPLE 

This section considers a simple example that illustrates several of 


the points discussed thus far. The potential caused by the charge 
density p is 


O(x.y,2) = a a ee a aay de, (1-10) 


where 


R=Y(x-2x'P +(y-y' P+ (2-7). (1-11) 


If the charge is confined to a thin flat plate in the plane z = 0 and 
has a surface charge density of p,, the potential in that plane is 


p(x’, y') 


cans (eet alae 


Consider now the special case where the plate is a square with sides 
a and the square is small enough that the charge density can be 
assumed to be constant over the square. Then 


(1-12) 


a/2 1 
o(x.y dx’ dy’, 
er rand a Bae EET y 
(1-13) 
and the potential at the center of the square is 
(1-14) 


af2 pa/2 1 
(0,0 dx’ dy’. 
a die = V(x" +(y')? . 
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_ 
x 


a/2 
FIGURE 1-7. Division of integration region into squares. 


e integral that we now want to evaluate numerically is 


a/2 e— 


af? Were fe + (yy 


we have normalized with respect to a so that INT is a pure 
ber. First let us take advantage of the two degrees of symmetry 
evaluate 


INT = Se hr dy'"=, (1-15) 


1 
0 V(t 


develop an algorithm for calculation of INT, divide the square 
‘into subsquares each with sides of h, as indicated in Fig. 1-7. The 
largest value of h that we can use is 


dxdy'. (1-16) 


a 
== 117 
n= a7) 
and in general 
a 
=— 1-18 
ve 2M° ( ) 
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where M is an integer. Over each subsquare we assume that we can 


teplace the integrand by its value at the center of that subsquare. 
Then (1-16) can be approximated by 


h2 


+ (1-19) 
miny 0.5h)° + (jh — 0.5h) 


M M 1 


2 
INT = : aa 
ud, % (i- 0.5) + (7-05) _ 


To complete the algorithm, we should utilize in succession values 
for M of 1,2,4,... up to some value determined by how much 
accuracy we wish in the final result and how much time we wish to 
spend in the computation. Also, we should apply Richardson ex- 
trapolation to the computed results. A simple computer program to 
carry out this process follows. 


FORTRAN Program for Evaluation of the Potential Integral 


program integral 

€ evaluation of potential integral 
integer M,i,j,k,N 
real $,Q(10,10) 


print *,'how many values of M?' 

read*,N 

print® 

print* Potential integral’ 

print*,' Evaluated with pulse functions' 
print* 
print*," “ 
print* 

M=1 

do k=1,N 


integral extrapolation’ 
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s=0 
do j=1,M 
do i=1," 
s=S+1/(sqrt(Ci-0.5)*Ci-0.5)+ 
(j-0.5)*(j-0.5)))/M*2. 
Qck,1)=S 
end do 
end do 
c=1 
do i=2,k 
C=Cx2 
Qk, i =CCHQ(K, 1-1) -ACK= 1,171) / 00-1) 
end do 
print*,M,@(k,1) ,Q(k,k) 
M=M*2 
end do 
end 


following table. 

M INT Extrapolations 

1 2.828427 2.828427 

2 3.150529 3.472631 

4 3.330882 3.524102 

8 3.426362 3.525561 

16 3.475469 3.525497 

32 3.500366 3.525495 

64 3.512901 3.525492 

128 3.519191 3.525496 


a value for M of 128, whereas the extrapolated values have 
arly converged to a value of 3.52549 by M of 16. The sixth 
cimal digit exhibits round-off error. This integral can be evalu- 
analytically and to six significant figures agrees with our 
ical evaluation. The integrand in the continuous integral has 
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a singularity at the origin, and the analytical solution must handle 
this in a limiting fashion. In our numerical solution, the points at 
which we evaluated the integrand did not include the origin and 
hence we avoided the problem. 


PROBLEMS 


1, Use of the Richardson extrapolation method with three com- 
puted results for h = 4s, 2s, and s leads to an extrapolation that 
can be expressed as 


E, = aC(s) + BC(2s) + yC(4s). 


What are the three coefficients a, B, and y? 


2. Suppose that a finite-difference computation for some quantity 
C has given the following results: 


h € 
1.0000 26.8285 
0.5000 16.2720 
0.2500 12.6726 
0.1250 11.1832 
0.0625 10.5062 


(a) What is your estimate of the true value of C? 
(b) How accurate do you think your value is? 


COMPUTER PROJECT 1-1 


Write and run a computer program to evaluate numerically the 
integrand 
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ITE DIFFERENCES IN ONE DIMENSION 


irst method we consider for solving differential equations is 
ly the simplest, the finite-difference method [1, 2]. The first 
n applying the finite-difference method is to select a discrete 
lues of x, the mesh points, as illustrated in Fig. 2-1 for one 
sion. If the spatial interval that separates the points in the 
is denoted by A, the mesh points can be denoted by the 
te index k, where 


ty = kh. (2-1) 


hen concentrates on the values of the function f at the mesh 
and in a sense ignores the values at other points. To simplify 
ition, One associates with the function f of the continuous 
x, the function F of the discrete index k, defined by 


P(K) = f(kh). (2-2) 


© next step in the solution is to replace the differential 
on, which involves derivatives of a function of the continuous 
€ x, by algebraic expressions involving the function of the 
te index k. First, we need approximations for the derivatives. 
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a oO is 2 


~ —+ — 


FIGURE 2-1. Mesh of points in one dimension. 
Consider the Taylor expansion 
he 
f(x +h) = L(x) + Af (x) + SF"(x) + net, (2:3) 


from which 


f(x +h) - h 
= : FO) pea) Ape, (2-4) 


The left side is an approximation to the derivative S'(x), but as h 
approaches 0, the error approaches 0 only as the first power of h 
and thus this is not a desirable approximation. The left side of (2-4) 
is a better approximation to the derivative at the midpoint. To see 
this, consider the expansions 


f(x +h) = 4x + 5] + : {x + 5] ca rls 5] 
h 
+B »(x + 3] - (2-5) 
and 
f(x) =s(r+ | = ‘ {x + 4] df r(x fs 5] 
P h 
~al"(s+5)+- (2-6) 
Subtraction of these expansions yields 
fie P ole 8 Male 8 on, 2m 
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that the error decreases as the second power of / and hence 
js a more desirable approximation. This is a good approxima- 
on for the first derivative between nodal points, but usually we 
d approximations at the nodal points. To find such an approxi- 
jon, combine the expansion 


h a 
aa = t) = f(x) — Af'(z) + SP") — Plz) + > (8) 
(2-3) to get 


f(a +h) — f(x —h) 


he 
7 =P (x) + ZP"a) tor. 29) 


is now the sort of approximation that we wish to use. In a 
ir Manner, One can show that 


a; = —h 1 
GHEE I= 27) P= pay 4 Petey gon, 


y (2-10) 


the left side of this equation is a good approximation for the 
nd derivative, In general we prefer to use approximations such 
(2-9) and (2-10), referred to as central difference formulas, 
because of the greater accuracy. 

Several generalizations are possible. The left side of (2-10) 
wolves what is known as a three-point formula, Even more accu- 
approximations involving five and more values are available, 
it it is questionable whether the greater accuracy is worth the 
sed programming time, and we shall not use them. Approxi- 
s for higher-order derivatives can be derived using the same 
oach, but approximations for first and second derivatives are 
it we need for the electromagnetic problems studied here. One 
eralization that we will encounter later in this chapter involves 
oximations for partial derivatives when we study problems in 
o and three dimensions. 
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2.2. ONE-DIMENSIONAL DIFFERENTIAL EQUATION EXAMPLE 
As a specific example, consider the equation 


af 

we tt4f=0 for0<x <1 (2-11) 
with the boundary conditions 

f(0)=0 and = f(4) = 1, (2-12) 


Suppose that we wish to use finite-difference approximations to find 
£(0.5). The largest value of h that we can use is 0.5. As shown in 
Fig. 2-2, there are three nodal values, of which only one, F()), is 
unknown. If we replace the derivative in (2-11) by the three-point 
central difference formula centered at k = 1, the differential equa- 
tion becomes 


F(0) ~ 2F(1) + F(2) 


+4F(1)=0 2-13 
(05) (1) (2-13) 

and the boundary conditions lead to 
F(0)=0 and F(2) =1, (2-14) 


Solution for F(1) gives 


F(0) + F(2) 


F(t) = 


(2-15) 


Insertion of the boundary values from (2-14) into (2-15) and evalua- 
tion of F(1) leads to 


F(1) = 1.0000. (2-16) 
i) 1 2 k 
to) 05 L # 


FIGURE 2-2. Three-point mesh for differential equation. 
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5 1 2 3 4k 
ri 05 LO & 


FIGURE 2-3. Five-point mesh for differential equation. 


is gives us a first value for f(0.5) of 1.0000. For a more accurate 
alue, take h as 


h = 0.25, (2-17) 
ind then we have the situation shown in Fig. 2-3, with two external 
points at which the functional values are determined by the 
dary conditions, namely 

F(0)=0 and F(4)=1 (2-18) 
d three internal points at which the functional values are to be 
termined. Use of the central-difference formula (2-10) at each of 
e internal mesh points leads to the three equations 


F(0) — 2F(1) + F(2) 


+4F(1)=0 
0.0625 
F(1) — 2F(2) + F(3) +4FQ)=0 (219) 
0.0625 
F(2) — 2F(3) + F(4) Aaah a 
0.0625 , 


Which can be rearranged as 


1.75F(1) — F(2) = F(0) 
— F(1) + 1.75F(2) — F(3) = 0 (2-20) 
— F(2) + 1.75F(3) = F(4). 


Substitution of the boundary conditions and direct solution of these 
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equations gives 


F(1) = 0.5378 
F(2) = 0.9412 (2-21) 
F(3) = 1.1092. 


The value for F(2) is an improved value for (0.5). We have 
calculated f(Q.5) for two values of A, 0.5 and 0.25. We can continue 
the sequence of calculations by using for h, the values 


Ee for M = 1,2,4, 2-22 
=r for M = 1,2,4,..., (2-22) 


where we have already used M = 1 and 2. The values of f(0.5) that 
result are shown in the following table. 


M £(0.5) 


1.0000 
0.9412 
0.9292 
0.9263 
0.9256 
0.9255 


NAeckRN 


l 
3 
Analytical solution of the differential equation gives 0.9254. 


2.3. FINITE DIFFERENCES IN TWO DIMENSIONS 


As we just illustrated in the one-dimensional example, the first step 
in the finite-difference method is to establish a mesh of points. 
Figure 2-4 illustrates a mesh in two dimensions. As in that illustra- 
tion, a uniform mesh is usually used with simple problems. Func- 
tions of discrete indices are then defined by 


F(i,j) = f(ih, jh). (2-23) 


If the complexity of the field varies over the region of interest, a 
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Gk+l G+Lkt) 
. . 


d-Lk) GK) i+1,k) 
. . . . . 
G@k-D 
. . . . . 


FIGURE 2-4. Two-dimensional mesh. 


iform mesh may be desirable. The derivative approximations 
oped for one dimension can easily be generalized to two and 
dimensions. For example, the finite-difference approximation 
Laplacian of a function of two variables is 


f(x +h, y)+f(x—-h,y) 
a Pf tf(x.y +h) + f(x,y —h) ~ 4f(x,y) 


ax? ej he : (28) 
in discrete terms is 
F(i— 1, J) + FL + 1,4) - 
V%i,/)- +F(i.j — 1) = +1) +4 F(t) (2:25) 
of this approximation with Laplace’s equation 
# ra 
= + a =0, (2-26) 


mple, results in the discrete equation 


—1,j)+F(i+1,j) +F(i,j-1) 
4 +F(i,j +1) —4F(i,j) =0. (2-27) 


Such equation is written for each unknown value of F, a set 
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FIGURE 2-5. Neumann boundary condition. 


of linear algebraic equations results. If the function f is specified 
on the boundary, a boundary condition known as the Dirichlet 
condition, the value of F is specified on all the boundary points. 
With this condition we wish to satisfy (2-27) at each internal point. 
If the internal point at which we wish to satisfy (2-27) adjoins a 
boundary point, some of the functional values in this equation are 
known. Another frequently encountered boundary condition is one 
known as the Neumann condition, where the normal derivative is 
specified and frequently to the value zero. For the situation of Fig. 
2-5, the point (i,j) is a boundary point where the Neumann 
condition 


(2-28) 


is to be satisfied. Because the functional value F(i, j) is unknown, 
we must use (2-27) at this point. This requires a value for F(i + 1, j). 
but the point (i + 1, j) is outside the computation region and thus 
the functional value F(i+1,j) is not available. The boundary 
condition (2-28) leads to 


F(i+1,j) =F(i- 1,7) (2.29) 


and we could use this relation to replace the value for the point 
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ie the computation region. Then (2-27) would be replaced by 
i — 1, /) + F(i, i + 1) + Fini — 1) — 4F(i,/) = 0. (230) 


‘The set of linear algebraic equations can be solved directly, but 
jike the one-dimensional example, even moderately small values 
h may result in so many unknowns that direct solution is not 
asible. Iterative, or indirect, methods have been found to be 
able for most two- and three-dimensional problems. For such 
ution, write (2-27) in the form 


F(i+1,j) + FU 1.i) + FGA + 1) + Fi, - 1) 
F(i,i) = 7 ‘ 


(231) 

The process of iteration starts with some initial values F “i, )) 
for the unknown values, and cycles repeatedly through all (i, i) 
that correspond to unknown values of F, and uses (2-31) in 
form 


Fi + 1,7) + Fi - 1,3) + Fi, 7 +1) + Fi, 7 — 1) 
i)= ; 


(2-32) 


an be used to calculate an improved set of values F\(i, j). Then a 
ond set of values can be calculated, and so on, using the general 


Fi + 1, jf) + FP" - 1,4) 
4F""\(i,j +1) + F"(i,7 - 1) 
a j 


Fei j)= (2-33) 
as before, the superscript indicates the number of the 
tion. This process can be continued until the solution has 
‘ged to suitable accuracy. Because the change from one step 
the iteration to the next may involve slight changes in the 
tional values, convergence is difficult to detect by examining 
cessive steps. A safer procedure is to double the number of 
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iterations successively and examine the changes after 1, 2, 4, 8, and 
So on, steps of iteration. When these numbers change little, it is 
safe to assume that the process has converged. 

The process of iterating b repeatedly using the same equations 
for the nodal values is known as relaxation, from carly structural 
applications in which the process was envisioned as “relaxing” 
stress. An improvement that speeds the convergence process is 
known as over relaxation, and for the example of Laplace’s equa- 
tion replaces (2-33) by 


F'(i,/) = F"\G,i) 
Fe“ +1, j) + F""\(i— 1,4) 
FPO MGI + D+F 'GI— 1) 


x, 
” 4 


— F*"(i,j) |. 


(2-34) 


If the constant R, known as the relaxation parameter, is equal to 
unity, the regular relaxation process results. Values of R greater 
than unity speed convergence, but too large a value of R results in 
numerical instability. For Laplace’s equation, a value of 2 gives 
instability, and values near 1.5 are shown to give good results. The 
program shown later uses a value of 1.4. 

In calculating the functional values for one iteration step, we 
used only the values from the preceding step, which is an iterative 
process known as the Jacobi process. For the first calculation in 
each step that is all we can do. For the second calculation, we have 
a choice, because the first calculation has given us a more accurate 
value for one of the numbers. The iterative process that always uses 
the most recent numbers available is known as the Gauss—Seidel 
method. The Gauss-Seidel method can be shown to converge faster 
for a wide range of problems, and in addition, is easier to program 
since we do not need to store more than one value for each 
unknown. 

The initial set of values F"(i, j) required to start the iterative 
process can be chosen arbitrarily as all zero or all unity. Conver- 
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to the desired solution occurs faster if one can start the 
with a set of initial values that are reasonably close to the 


values. 


4. A TWO-DIMENSIONAL CAPACITANCE EXAMPLE 


mple two-dimensional example is the quasistatic analysis of a 
d parallel plate transmission line, the cross section of which 
in Fig. 2-6. The potential satisfies Laplace’s equation 


2, 2, 
Ee 2S ny (2-35) 
ax? ay? 
convenience, we will take the potentials on the upper and lower 
as +1 and —1, respectively, and the potential on the outer 
d as 0. To minimize the computational effort, we can take 
ntage of the top-to-bottom antisymmetry and the left-to-right 
netry to replace the situation by the one shown in Fig. 2-7. The 


——$<$<<= 


SURE 2-6. Cross section of shielded parallel plate transmission line. 


FIGURE 2-7. One-fourth of cross section. 
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methods discussed in the preceding section can be used to solve for 
the potential (x, y) and then in turn we can calculate the capaci- 
tance. The fundamental steps in a computer algorithm to carry out 
this computation are: 


j. Establish parameter values. 

2. Establish boundary conditions and initial values. 
3. Iterate to determine the potential values. 

4. Calculate the capacitance. 


We need first to select values for the parameters a, b, w, and d. 
Since the capacitance is not affected by scaling, we will use instead 
the three normalized parameters a/d, b/d, and w/d. Then we 
establish the potential on the shield of 0 volts and on the center 
conductor of 1 volt. Next we perform a set of iterations. To perform 
one iteration, we raster scan through all points at which the 
potential is unknown and at each point utilize the overrelaxation 
expression. As discussed in the preceding section, a good approach 
is to examine the results after one iteration, two iterations, four 
iterations, and so on. After each such set of iterations, we calculate 
the capacitance and compare it to the previous value. We stop the 
iteration after the absolute difference of the two capacitance values 
is less than some tolerance value that we select as one of the input 
parameters, 

For each set of potential values, we calculate the capacitance. 
The capacitance can be expressed in terms of the charge on the 
center conductor as 


Q 
c=s, (2-36) 


where the voltage V is the difference between the potential values 
on the two conductors. Because this potential difference is 2 volts, 


2 
2 volts © 


(2-37) 


The charge QO can be calculated from the two-dimensional version 
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f Gauss’s law as 

O = eq f Eau dl, (2-38) 
the line integral surrounds the center conductor and the 


ubscript on E means the normal component away from the center 
ductor. Because 


E. 


‘oul 


ab 


integral can be computed as 


0=«f (3). dl, (2-40) 
0 ab 
c= sot (Fe), dl. (2-41) 


ince we are analyzing just part of the total region, we will evaluate 
jis line integral for that subregion and multiply the result by 2. 
_ To complete the algorithm, we need to select values of the 


finite-difference distance A and select a criterion for stopping the 
eration. The largest value that we can take for h is 


ha = 5 (2-42) 


Jas te zy (2-43) 
M 2M 

vhere M is an integer whose value we shall, in succession, take as 

:2,4,... up to as many values as we wish to use. 


32° FINITE-DIFFERENCE METHOD 


We thus calculate a value of capacitance for each value of h. The 
final step is to apply Richardson extrapolation, as discussed in 
Chapter 1, to these capacitance values to arrive at our best esti- 
mate. A FORTRAN program that implements the complete algo- 


rithm follows. 


FORTRAN Program for the Capacitance of a Shielded 
Parallel-Plate Transmission Line 


program shcapfd 
capacitance of shielded capacitor 
calculated with finite differences 
Line integral of E field 
solution by iteration 
real add,bdd,R, delta,cap,old 
real z 
integer C,W,WW,M,nx,ny,i,j,k,kk 
real V(200,200),Z2(10,10),¥(10) 
C input parameters 
print* 
print* 
print* 
print* 
print*,' shielded capacitor’ 
print*,' finite-difference approximation’ 
print*,' capacitance from Line integral’ 
print*,' solution by iteration’ 
print* 
print*,' what is a/d?' 
read*,add 
print* 
print*,'what is b/d?! 
read*,bdd 
print* 
print*,'"how many values of M?' 
read*,WW 
print* 


a0a00 
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print*,' enter tolerance’ 
read*,tol 

print*® 

print® 

R=1.4 

print parameters 

print*,' shielded capacitor’ 
print*,' a/d is ',add 
print*,' b/d is ',bdd 
print*,' tolerance is ',tol 
print* 

print*,' M original cap extrapolations' 
print* 

start calculation 

M=1 

w=1 

do while (W.LE.WW) 

nx=add*M 


_ ny=bdd*™ 


set initial and boundary conditions 
doj=0,ny 
doi=0,nx 
VCi,j)=0 
end do 
end do 
do j=1,ny-1 
do i=0,nx-1 
VCi,j)=1 
end do 
end do 
cap=1 
old=0 
k=1 
kk=1 


do while (abs(old-cap).G6T.tol) 


old=cap 
do while Ck.LE-kk) 
doj=1,ny-1 
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doi=0,nx-1 
if(j -EQ@. M) then 
if Ci -GT. M) then 
Zz=VCi-1,j)+VCi41,j5) 
z=z4+VCi,j-1)4VCi,j+1) 
delta=z/4-V(i,j) 
VCi,j)=VCi,j)+Redelta 
end if 
else 
if (i.E@.0) then 
z=24VCi+1,j) +VCi,j-1) 
z=z+VCi,j+1) 
delta=z/4-vCi,j) 
VCi,j=VCi,j) +R¥delta 
else 
z=VCi-1,j)4+VCi4+1,j5) 
z=z4+VCi,j-1)+VCi,j+1) 
delta=z/4-VCi,j) 
VCi,j)=VCi,j)+Redelta 
end if 
end if 
end do 
end do 
k=k+1 
end do 
2z=(VC0,1)+VCO,ny-1))/2 
do 7=1,nx-1 
Z=24+VCi,1)4VCi,ny-1) 


end do 
do j=1,ny-1 
z=z+V(nx-1,3) 
end do 
cap=8.854187*z 
kk=2*kk 
end do 


ZZ(W,1)=cap 


€ extrapolation calculation 


net 
do j=2,W 
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c=2*C 

ZZCj ,W)=CC*ZZ0 5-1, W)-ZZCj-1,W- 129 / 0C-1) 
end do 

print results 

print*,YCW),ZZ(1,W) ,2Z(W,W) 

W=W+1 


results of a computation run for the parameter values 


a) 
= 
b 
—= 2-44 
==2 (2-44) 
“ 1 
= 
ire shown in the following table. 

M Capacitance Extrapolations 

1 35.41675 35.41675 

2 29.93558 24.45442 

4 27.83906 26.17192 

8 26.89887 26.01054 


OPEN REGIONS 


. use of the finite-difference method to solve problems such 
€ example in the preceding section is straightforward. The fact 
we need to consider only the region inside the shield means 
we have only a finite number of mesh points. Situations with 
regions, however, present a serious problem, because there 
an infinite number of mesh points, which clearly cannot be 
idled. One approach to such a situation is to restrict the compu- 
tion to a finite portion of the space and then to successively 
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expand this portion of the space. For example, if we consider the 
parallel-plate transmission line, we can retain the shield but make 
the shielded region larger and larger. Richardson extrapolation can 
then be utilized to estimate the value of capacitance with the shield 
at an infinite distance. The following table gives the results of such 
a computation. 


a/d =b/d Capacitance (pF /m) Extrapolations 
2 26.01054 26.01054 
4 20.18555 14.36056 
8 19.07868 19.17556 
16 18.81944 18.69644 


Although it is not clear from this table how well the value has 
converged, the approximate value of 18.7 does agree with the value 
we will derive in a later chapter from the much more efficient 
moment method. Although this method of handling open regions 
does work, it leads quickly to very large numbers of points and is 
not very computationally efficient. The moment method that we will 
treat later proves to be much better suited to open regions. 


2.6. GENERALIZATIONS 


The replacement of derivatives by finite-difference approximations 
clearly can be applied to other partial differential equations. Al- 
though we here discussed only one- and two-dimensional problems, 
the generalization to three dimensions presents no problems in 
principle. The additional dimension does, however, lead to a larger 
number of points and increased run time. This is why the tendency 
is to solve two-dimensional problems wherever possible. Three- 
dimensional problems will be discussed further in Chapter 4, where 
time-domain problems are discussed. 

Because the finite-difference method is based on a rectangular 
mesh, rectangular conditions on rectangular boundaries can be 
satisfied in a natural way. When the boundary is curved, special 
modifications must be made. One way is to use straight “stair-step” 
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dary approximation inside or outside the true boundary, as 
n in Fig. 2-8. One can solve the problem twice, once with an 
approximation and once with an outside approximation, and 
average the results. A more sophisticated way to accommodate 
d boundaries is to modify the finite-difference approximation 
derivatives to allow unequal spacing whenever a nearest point 
on the boundary, as shown in Fig. 2-9. Although there may be 


FIGURE 2-9. Unequal spacing for curved boundary nodes. 
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problems where the additional programming effort is warranted, 
usually a stair-step approximation is adequate. 


PROBLEMS 


1, The simple one-dimensional ordinary differential equation 


df 
ae tf? 


with the boundary conditions 


f(0)=0 and f(i)=1 


is to be solved with the finite-difference method. For a very 
crude mesh, take A = 0.5. What is (0.5)? 


2. For the values of at the nonuniformly spaced points in the 
figure, calculate an approximate value for 4° /dx? + 47 /ay?. 


COMPUTER PROJECT 2-1 


This computer project is the calculation of the capacitance per unit 
length for the square coaxial transmission line shown in the figure. 
Use the finite-difference method and calculate the value to an 
accuracy of 0.1 pF/m. Your report should include: 


a, A statement of the problem. 

b. Equations and algorithms used. 
c. A program listing. 

d. Intermediate numerical results. 
e. The answer. 
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PUTER PROJECT 2-1. Cross section of sqaure coaxial-like transmis- 
line. 
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inite-Difference 
termination of Eigenvalues 


EIGENVALUES IN ONE DIMENSION 


extremely important problem in electromagnetics is the deter- 
ation of resonant frequencies of a structure such as a cavity. 
thematically, this is described as an eigenvalue problem. A sim- 
eigenvalue problem in one dimension involves determining the 
es of & that are consistent with nontrivial solutions of the 
erential equation for f(x): 


—5+kf=0 with f(0)=0 and f(1)=0. (+1) 


general solution of the equation is 
f(x) =A cos kx + Bsin kx, (3-2) 


id the imposition of the boundary conditions requires A = 0 and 
na, where n is an integer. One way to find the values of k 
nerically is to use the finite-difference approximation for the 
ond derivative, which leads to the set of equations 


fin — 2fi + fiat 


= + kf, =0, (3-3) 
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which can be written in the form 


fiat ef, fir = 0, (3-4) 


where 
@=2—h?k?. (3-5) 


Nontrivial solutions of these algebraic equations exist only if the 
determinant of the system of equations is zero. When the number 
of unknown functional values is small, the determinant can be 
evaluated explicitly. For example, if 4 = 0.25, there are three 
unknowns, and setting the determinant equal to zero gives 


@ 1 0 
=0. (3-6) 


Evaluation of the determinant leads to the equation 
a(a? — 2) =0, (3-7) 
which has the solutions 
a@=0,4+V2 (3-8) 
which leads to k values of 
k = 3.0615, 5.6569, and 7.3910. (3-9) 
Compared to the true solutions of 3.1416, 6.2832, and 9.4248, the 
first value is a fairly good approximation, the second is so-so, and 
the third is not a very good approximation. 
3.2, WAVEGUIDE-MODE EXAMPLE 
A more realistic example is the problem of finding in two or three 


dimensions the eigenvalues that determine the propagation param- 
eters of waveguide modes, resonant frequencies of resonators, and 
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i i le, consider a 
other physical parameters. As a simple example, 
gular Rerncaibie, with propagation in the z-direction. For the 
no the four field components E,, E,. Ay, and H, can be 
in terms of E,. In turn, E, can be written as 


E,(x, y,Z,t) = E(x, yeh), (3-10) 
e E satisfies 
PE PE ’ 
gee ye tea (3-11) 
E=0 (3-12) 


boundary. The parameter k determines the phase parameter 
ough 


k? = w*pe — B?. (3-13) 


eigenvalue problem is to solve for E(x, y) and k simultane- 

so that (3-11) is satisfied. - 

use the finite-difference method, a mesh of points is estab- 
separated by the uniform interval of A. With the previously 
approximation for the Laplacian, the discrete equation that 


he 
+k°E(i, j) = 0, (3-14) 


(i + 1, j) — E(i—1,j) — E(i,f +1) — Eli.f —1) 
+(4—K7n)E(i, 7) =0. (3-15) 
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With the notation 
a=4—k*h?, (3-16) 
this equation is 
~E(i + 1,7) — Ei - 1.) — E(i,j + 1) - E(i,j - 1) 
+a@E(i,j)=0. (3-17) 


This is a set of linear algebraic equations, one for each interior 
point, These equations can be written in matrix form as 


A(a)E =0, (3-18) 

where the dependence of the matrix coefficients on @ i ici 
: 2 is explicitl 
written. For the set of equations (3-18) to have nonzero foliitene 
for E, the determinant of the matrix must satisfy 
det[ A(a)] = 0. (3-19) 


As is well known, this determinant is a polynomial in a with degree 
equal to the number of points. For 


(3-20) 


FIGURE 3-1. Mesh of points for finite-differences with h = a/4. 
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equations (3-18) are 


aE, —E,=0 
—E,+aE,—E,=0 
E,+ aE, =0. (3-21) 


determinant of this set of equations, which happens to be the 
determinant as that considered in the preceding section, 


a -l 0 
—1 a —1|)=0, (3-22) 

0 —1 a 

ch leads to the polynomial equation 
a —2a=0. (3-23) 


polynomial has the zeros — V2, 0, and ¥2, which lead to 

for ka of 6.4322, 8.0000, and 9.3074. These are approxima- 
to the first three eigenvalues. To get more accurate values, we 
of course to use smaller values of A. A value for h of 


A= 3-24) 
=— (3- 


example, leads to the situation shown in Fig. 3-2, which shows 
there are 21 internal points. To write explicitly the determinant 
corresponds to (3-22) as a polynomial in « is not an attractive 


FIGURE 3-2. Mesh of points for finite-differences with h = a/8. 
? 
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task and for even smaller values of 4 almost impossible. A number 
of algebraic techniques [1] have been devised to find eigenvalues of 
such determinants. In the next two sections we consider two numer- 
ical techniques. 


nt of a diagonal matrix is the product of the diagonal 
the determinant of U is unity. Thus, with this algorithm, 


det[ A] = det[L], (3-31) 


since L also is a diagonal matrix, its determinant is the product 
diagonal terms, and thus 


det[ A] = TH. i). (3-32) 
- 


3.3. NUMERICAL EVALUATION OF THE DETERMINANT 


One approach is to use an algorithm that numerically evaluates the 
determinant [2]. First, we need to count and number the internal 
points. Then we need to evaluate the matrix coefficients that 
correspond to (3-19), Evaluation of the determinant can be carried 
out using a modification of an algorithm for matrix inversion known 
as LU factorization. If the matrix A can be factored as 


A=LU, (3-25) 


Although unlikely, one of the values of ka chosen could be a 
of the determinant (and in the first example we will soon 
ider, this happens). Provision must be included to recognize 
possibility and avoid trying to divide by zero in the factorization 


cess. A FORTRAN program for the TM modes of a rectangular 


where L is a lower diagonal matrix with zeros for all elements 
above the diagonal and U is an upper matrix with zeros for all 
elements below the diagonal. The principal use of this approach is 
the solution of the matrix equation 


ide follows. 


program tmdet 
integer M,R,RR(100,100) ,W,F 


integer i,j,k,N,imax,indx(110) 


can (3-26) real vv¢110) 
real aamax, sum, dum, D 
deyrteltiny/ thts input parameters 
print*,'what is M?' 
pane G27) read*,M 
i j ini Q" 
and then writing this as the two equations prints, ‘what is minimum ka? 
read*,minka 
LY=B (3-28) print*,'what is maximum ka?! 
and read*,maxka 
UX =Y. (3-29) count and number internal points 


We wish to use this approach to evaluate the determinant of the 


R=0 
doj=1,2*M-1 


matrix A. The determinants that correspond to these matrices are ol Sal adel 
related by Bae 
RRCi,jI=R 
det[ A] = det{ L]det[U]. (3-30) end do 
end do 
One factorization algorithm, known as Crout’s algorithm, results in N=R 


the matrix U having ones for the diagonal terms and since the 


Print*,'number of internal points 


real alpha,ka,kh,A(110,110) ,minka,maxka 


is',N 
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print column headings 


print* 

print*,' ka sign(det)! 
print* 

calculate values of ka and alpha 
do W=0,10 


ka=minkatW*(maxka-minka)/10. 
kh=ka/4./M 
alpha=4.-kh*kh 
do j=1,N 
do i=1,N 
ACi,j)=0. 
end do 
end do 
evaluate matrix elements 
do j=1,2*M~-1 
do 1=1,4*M-1 
ACRRCI,j),RRCi,J))=alpha 
ACRRCi,j),RRCi-1,j))=-1. 
ACRRCI,j),RRCI+1,j))=-1. 
ACRRCi,j),RRCI,-J-1))=-1. 
ACRRCI,j),RRCi,j+1))=-1. 
if (j.E@.1) then 
ACRRCi,5),RRCi,d-1))=0. 
else if (j.£Q.2*M-1) then 
ACRRCi,j),RRCi,j+1))=0. 
end if 
if Ci-E@.1) then 
ACRRCi,j),RRCI-1,j))=-1. 
else if (i.EQ.4*M-1) then 
ACRRCi,j) ,RRCI+1,j))=-1. 
end if 
end do 
end do 
start LU decomposition 
D=1. 
Fe1 
do 7=1,N 
vu Cid=t. 
end do 
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do i=1,N 
aamax=0. 
do j=1,N 
if Cabs(ACi,j)).GT-aamax) aamax= 
abs(ACi,j)) 
end do 
if Caamax.E@.0.) then 
F=0 
else 
vv (i)=1./aamax 
end if 
end do 
do j=1,N 
do i=1,j-1 
sum=ACi,j) 
do k=1,i-1 
sum=sum-ACi,k)*ACkK,j) 
end do 
ACi,j)=sum 
end do 
aamax=0 
do i=j,N 
sum=ACi,j) 
do k=1,j-1 
sum=sum-ACi,k)*ACkK,j) 
end do 
ACi,j)=sum 
dum=vv(i)*abs (sum) 
if (vv(i).GT.aamax) then 


imax=i 
aamax=dum 
end if 
end do 
if (j.NE.imax) then 
do k=1,N 
dum=ACimax,k) 
ACimax,k)=ACj,k) 
ACj,k)=dum 
end do 
D=-D 
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with large matrices. Then one repeats the process for another 

s, spanning the region where the determinant changed sign 
first calculation. In this way, a new decimal digit is deter- 
d, and the process can be continued until the accuracy desired 
‘ed. As an example, for M = 2, the following results are 


vv Cimaxd=vv¢j) 
end if 
indx(j)=imax 
if (ACj,j).£@.0) then 


F=0 
else ed. 
if (j.NE.N) then “ 
dum=1./ACj,j) ka sign(det) 
do i=j+1,N 0 = 
ACi,j)=ACi,j)*dum 1 =k 
end do Z =T 
end if 3 -1 
end if 4 =f 
end do 6 a 
c calculate sign of determinant 7 #E 
if (F.E@.0) then 8 0 
D=0. 9 =1 
else 10 1 
do j=1,N 
D=D*ACJ,j)/abs(A(j,j)) ‘can see that there is a zero between 6 and 7, one at 8 exactly, 
ae one between 9 and 10. If we desire to get a more accurate 
end if of the first one, for example, we compute ka between 6 and 7 
c print solution for one value of ka get the following results. 
print *,ka,D 
c repeat calculation for another value of ka ka sign(det) 
end do 
ais ela 
An algorithm that evaluates the determinant can be used in “ =. 
several ways. The determinant is a function of a, or equivalently of 6. 4 -1 
ka, and the Newton—Raphson method for finding zeros of functions 6: 5 ot 
can be combined with the determinant evaluation algorithm. The 66 44 
bisection method is another useful way to find zeros of functions. 67 +1 
For interactive computations, a scheme based on decimal division 68 #1 
seem to more useful. The determinant is evaluated at 11 values of 6.9 41 
ka, say from 0) to 10, and one looks for successive values of ka for 70 ay 


which the sign of the determinant changes sign. If this process is to 
be followed, only the sign of the determinant must be calculated 


and not the value. This simplifies numerical problems that may zero is seen to be between 6.4 and 6.5. After several more such 
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calculations, we get a value of approximately 6.43215. For three 
values of M, the first three zeros are given in the following table. 
The table also shows the extrapolated value and the analytic result. 


M 
1 Zz 4 Extrapolation Analytic Result 
6.43215 6.87265 6.98655 7.02955 7.02481 
8.00000 8.65915 8.82875 8.89168 8.88577 
9.30735 10.79385  11.19205 11.36022 11.32717 


This entire process could be programmed to completely automate 
calculation of the eigenvalues. Although in principle this approach 
should work as well for the TE modes as it does for the TM modes, 
unfortunately some of the zeros of the determinant are multiple 
and are extremely difficult to locate. 

A number of other matrix-based algorithms have been developed 
for calculation of eigenvalues. The well-known power method is an 
iterative process that converges to the largest eigenvalue. To deter- 
mine the smallest eigenvalue, which is what we usually want, the 
process can be applied to the inverse matrix. A more elaborate 
algorithm that is harder to understand is the QR method, which in 
principle produces all the cigenvalues. The reader is referred to the 
literature [2] on algebraic algorithms for details. 


3.4, ITERATIVE SOLUTION METHODS 


General algebraic methods and methods such as that in the preced- 
ing section exist for finding the eigenvalues by direct calculation 
from the matrix representation. However, we clearly need a better 
method for very large matrices, and iterative methods are one such 
approach. To utilize an iterative method, we need equations to 
determine E(i,j) and k. The equation for E(i, j) can be derived 
from (3.9) as 


_ E(i+1,j) + E(i—1,j) + Eli,i + 1) + Eli,i — 1) 


E(i,j) a 


(3-33) 
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uld 
t zero. If the values of E were known exactly, (3-15) cor 
ees if a is not zero. If the values of E were known exactly, 
could be solved for k* to give 
WE 
= 


= (3-34) 


value of (i,j). A useful equation which, in a sense, is a 
d average can be derived by multiplication of (3-11) by E 
integration over the entire region to give 


[[ BV?Ededy +k? ff EP dedy =0. (3-35) 
Solution for k? yields 
2, 
pn hacks ——. (3-36) 


i 
of the finite-difference approximation in this expression gives 


¥Y LEG, NE + 1,7) + EG - 1,4) + EGF + 1) 
ae: +E(i,j - 1) — 4B, J) 


‘Sg LE lei I 
(3-37) 
of (3.16) simplifies this equation to 
LY LEW MEG + 1,4) + BG - 1) 
_ : +E(i,f+1)+E(,i- DI G38) 


LLG. 


Alternate use of (3.33) and (3.38) should converge to a solution 
E and k. Which mode is converged to depends on the initial 
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conditions. For example, with the initial condition that E is set to 
unity at all interior points, the convergence is usually to the lowest 
mode. Many initial conditions seem to converge to the lowest mode, 
To converge to other modes, one must select carefully the initia] 
conditions, and this is probably not practical. A FORTRAN pro- 
gram for TM modes that uses unity for initial conditions follows. 


program wgotm 

Cc iterative solution for eigenvalue of TM 
mode for rectangular waveguide 
real (130,130) ,Q(10,10) 
real R,tol,num,den,beta,ka,oldka 
integer M,i,j,n,nn,W,WW,C,MXx 
print*,'how many values of M?' 
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num=0 
den=0 
do j=1,2*M-1 
do i=1,4*N-1 
beta=f(it+1,j)+f(i-1,j) 
beta=beta+f(i,j+1)+f(i,j-1) 
num-num+f(i,j)*beta 
den=dentf(i,j)*fCi,j) 
end do 
end do 
alpha=num/den 
ka=4*M*sqrt(4-alpha) 
do j=1,2*M-1 
do 1=1,4*M-1 


read* ,WW beta=fCit1,j)+f(i-1,j) 
print*,'what is tolerance?’ beta=betatf(i,j+1)+f(i,j-1) 
read*,tol Ci, j=FCi,j)+R*(beta/alpha- 
R=0.7 #(i,j)) 
M=1 end do 
do W=1,WW end do 
print*,' iterations alpha kat n=nt1 
print* end do 
doj=0,2*™ print*,(n-1),alpha,ka 
doi=0,4*mM nn=2*nn 
#4, j)=0. end do 
end do Q(1,W)=ka 
end do c=1 
do j=1,2*M-1 do i=2,W 
do i=1,4*M-1 c=24*C 
fCi,j)=1. QCi,W)=CC*QC5-1,W)-QCI-1,W-199/ 00-1) 
end do end do 
end do do 71=1,25 
n=1 print*® 
nn=1 end do 
| ka=1 print*,' TM mode of rectangular 
oldka=0 waveguide" 
| do while (abs(ka-oldka).GT.tol) print*,' jterative finite-difference 
oldka=ka solution' 
do while (n.LE-nn) print*,* tolerance is',tol 
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print® 
print*,' mM ka extrapolation’ 
print* 
Mx=1 
do j=1,W 
print*,Mx,@(1,j),Q¢j,j) 
Mx=2*Mx 
end do 
M=2*M 
end do 
end 


PROBLEMS 
1. Use the finite-difference method to find an approximate value of 


the lowest normalized eigenvalue ka for the TM mode defined 


by 
VE, + KE, = 0, 


with E. = 0 on the boundary, for the mesh shown in the figure. 
Compare with the analytic solution. 


o 
i) 
wir 
o 


2 
PROBLEM 3-1. Mesh for eigenvalue calculation. 


2. Use the finite-difference method to find an approximate value of 
the lowest normalized eigenvalue ka for the waveguide TM 
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VE. + k?E, = 0, 


th E. = 0 on the boundary, for the mesh shown in the figure. 
that the spatial intervals h, and h, are not equal. Compare 
th the analytic solution. 


PROBLEM 3-2. Mesh for cigenvalue calculation. 


the finite-difference method to find approximate values of 
normalized eigenvalue ka for the nonrectangular waveguide 


PROBLEM 3-3. Mesh for eigenvalue calculation. 
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TM mode defined by 
VE, + k*E, = 0, 


with E. = 0 on the boundary, for the mesh shown in the figure. 
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ite-Difference 
e-Domain Method 


|. WAVE EQUATION IN ONE SPATIAL DIMENSION 


e finite-difference time-domain method is especially useful for 
lution of wave equations. Let us begin with the simplest 
n, that of one spatial dimension, where the wave equation is 


OE 
ar 


PE 


ees 4-1 

ax~ eH 
general solution of this equation consists of two waves, one 
g in the positive x-direction, so 


E(x,t)=E,(x — ct) + E(x + ct). (4-2) 
o functions on the right side of (4-2) are determined by the 
and initial conditions. The boundary conditions require 
edge of 


E(0.t)=b(t) and = E(a,r) =e(r) (4-3) 
¢ initial conditions are 

a 
E(x,0)=f(x) — and ap Els0) = a(x). (4-4) 
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If the finite-difference approach is used with the spatial variable x, 
and we take h as 


a 
h= NV’ (4-5) 
the central difference for the second derivative leads to 
: h? d?E(i,t) 
B+ 1,t) — 2E(i, 4) + E(t - 1,0) = a (4-6) 


dt 


where the partial differential equation now has been replaced by a 
set of ordinary differential equations. There are a total of N +1 
points, at two of which E£ is known from the boundary conditions, 
and N— 1 locations at which E must be determined from the 
N — 1 ordinary differential equations defined by (4-6). Although it 
is sometimes appropriate to solve these ordinary differential equa- 
tions analytically, usually a numerical solution is more useful. 


4,2. TIME QUANTIZATION 


In addition to quantizing the spatial variables, the time variable can 
be quantized in the same way using a time interval 5. The field 
variable as a function of discrete indices we now denote by E(i, n), 
and application of the central difference approximation to the 
derivative with respect to time leads to 


E(i,n + 1) — 2E(i,n) + E(i,n ~ 1) 


ee 
eh dnt) — BBE A) + EG — Ln) (47) 
ie 
Some rearrangement gives 
E(i.n +1) — 2E(i,n) + E(i.n — 1) 
=p [E(i+ 1,n) —2F(i,n) +F(i-1,n)], (48) 
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pay (4-9) 


onvenient to rearrange this equation for computation into the 


i,m + 1) = (2 — 2p?)E(i.n) — E(i,n — 1) + p?S(i,n), (4-10) 


S(i,n) = E(i-1,n) + E(i+ 1,n). (4-11) 


determined at the endpoints from the given boundary condi- 


E(0,n) = b(nd) and E(N,n) =c(nd). (4-12) 
> E is determined for n = 0 and for » = 1 from the initial 
itions, later values can be calculated explicitly from (4-10) and 
. It can be shown that this numerical process is stable if 
p<i. (4-13) 
we wish to investigate the conditions under which the 

rence approximations that we have used would be ex- 
d to be reasonably accurate. To investigate these questions, let 
nsider sinusoidal solutions. One set of continuous solutions 
the form 


o 
E(x,t)= sin( = ot}. (4-14) 
discrete solution can be described by 
o 
E(i,n) = sin( “ih = nad). (4-15) 
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Let us investigate the circumstances under which the finite-dif- 
ference approximation for the second derivative with respect to x js 
valid. From the continuous expression (4-14), 


PE o\? [@ 
fil alta, 


: (4-16) 


and for the discrete solution (4-15), some trigonometric manipula- 
tion leads to 


E(i + 1,n) — 2E(i,n) + E(i- 1,n) 
he 


2 oh _ [ox 
- (co = 1sin( “* - at}. 


The approximation for the second derivative clearly is valid if the 
approximation 


(4-17) 


wh , wh? 
ak ~~ Det 


(4-18) 


is valid. This approximation is good to 1 percent if wh/e is limited 
to approximately 


oh 
(=) = 0.35. (4-19) 
max 
In terms of the wavelength A, this limit is approximately 
h 
i* 0.055, (4-20) 


which means that we should have at least 18 points per wavelength, 
which seems reasonable. A similar calculation for the approxima- 
tion to the time derivative leads to 


(@8) max = 0.35. (4-21) 
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@8 = pkh (4-22) 


must be less than 1 for stability, (4-21) is satisfied automati- 

when (4-19) is satisfied. The conditions just derived show when 

would expect a reasonably accurate representation of the con- 
is solution. 


INITIAL CONDITIONS 


addition to the differential equation, two boundary conditions 
for all time ¢ and two initial conditions defined for all x are 
d to specify the solution. The two initial conditions are 
ly taken as the value of the function for time zero and the 


function into conditions on the discrete function is given by 
and (4-4). Use of the first derivative initial condition, however, 


aE E(x, 8) — E(x,0) 
oo es (4-23) 


but as we saw earlier, this approximation is accurate only to the first 
der in 6. A second-order approximation can be found from the 
or series expansion 


dE 2 a2 
E(x,5) = E(x,0) + 35 -(0) + Sar (440) +. (4-24) 
titution from the wave equation gives 
252 a, 


= dE cé 
E(x,5) = E(x,0) + 85 -(2,0) + “a ania (4-25) 


troducing the finite-difference approximation for the derivative 
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with respect to x yields 


= 
E(x,8) = E(x,0) + = = (2, 0) + oF 


x [E(x - el - pen ee +h,0)] +---, 
(4-26) 


which can be rearranged to give 
2 
E(x, 8) = (1 — p?)E(x,0) + 8¢(x) + (2,0), (4-27) 


where 


S(x,0) = E(x —h,0) + E(x +h, 0). (4-28) 


4.4. WAVES IN TWO AND THREE SPATIAL DIMENSIONS 


Although we have concentrated so far on the situation with one 
spatial dimension, the same equations with appropriate modifica- 
tions hold for two and three spatial dimensions. The finite-dif- 
ference approximations for the derivatives with respect to y and z 
follow the same pattern. The relation for one spatial dimension, 
(4-12), becomes 


F(i,j,n + 1) = (2 — 4p?) E(i, j,n) 
— E(i,j,n —1) + p?S(i,j,n), (4-29) 
where 
S(i,i,n) = E(i- 1, i,n) + E(i + 1,i,n) 
+E(i,j—1jn)+E(i,i+1,n) (4-30) 


and 


e=>- (4-31) 
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e condition in two dimensions for stability is 


psd, (4-32) 


is approximately 


p < 0.707. (4-33) 
initial condition becomes 


E(x, y,0) = f(x,y), (4-34) 


and the equation for E(x, y, 5) following the same pattern as the 


dimensional case becomes 


2 
E(x, ¥,8) = (I~ 2p?) E(x, ¥.0) + a(x, y) + F5(xy,0). 
(4-35) 


As an example, solution for the TM waves in a rectangular 
guide leads to a program like that shown below (which does 
include a provision for graphical output). 


Program for Calculating TM Wave in Rectangular 
ide Using Finite-Difference Time-Domain Method 
for Calculating the Frequency Spectrum 


Program wgtm 

calculates two-dimensional TM wave and its 
frequency spectrum 

for rectangular waveguide 

integer i,j,k,M,n,nt,nf 

real E(0:64,0:32),F(0:64,0:32) 

real 6(0:64,0:32) 

real rho,z,R(0:3200) 

real alpha,d,ad,Re,Im,Mag(1600) ,pi ,w df 


FINITE-DIFFERENCE TIME-DOMAIN METHOD 


pi=4.* atan(1.0) 
rho=sqrt(0.5) 
input parameters 
print*,'what is M??* 
read *,M 
print*,'what is initial distribution param- 
eter alpha?" 
read*,alpha 
print*,'how many time points?' 
read*,nt 
print*,'how many frequency points?! 
read*,nf 
df=4.0*M/nt/rho 
establish boundary conditions 
do i=0,4*M 
do j=0,2*" 
GCi,j)=0. 
FCi,j 
ECi,j)=0. 
end do 
end do 
establish initial conditions 
do j=1,4*M=1 
do j=1,2*M-1 
d=sqrt(Ci-M)*Ci-M)+(Cj-M)*Cj-M)) 
ad=alpha*d 
if (ad.LT.5) then 
GCi,j)=exp(-ad*ad) 
else 
G(i,j)=0. 
end if 
end do 
end do 
RCO)=G(3*M,M) 
calculate values for n=1 
do i=1,4 *M-1 
do j=1,2 * M-1 
S=GCi-1,j) +6 Ci+1,j)+6C1,j-1)4+6C7,j+1) 
z=(1-2*rho* ho)*GCi,i) 
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z=z+rho*rho/2.* $ 
FCi,j=z 
end do 
end do 
RC1)=FC3*M,M) 
print first two values of response 
do n=0,1 
print*,n,R(n) 
end do 
calculate for remaining values of n 
do n=2,nt~1 
do i=1,4*M-1 
do j=1,2 *M-1 
S=FCi-1,j) +FCi+1,j)+FCi,j-1)+ 
FCi,j+1) 
z=(2-4*rho*erho)* FCi,j)-GCi,j) 
z=z+rho*rho*e S 
ECi,j)=z 
end do 
end do 
R¢n)=EC3*M,M) 
print*,n,R(n) 
doi=1,4*M-1 
doj=1,2*M-1 
GCi,j=FG,j) 
FCi,j=ECi,j) 
end do 
end do 
end do 
calculate frequency spectrum 
do k=0,nf-1 
Re=0. 
Im=0. 
w=2.*pi/nt * k 
do n=0,nt-1 
Re=RetR(n)*cos (wen) 
Im=Im+#R(n)*sin(wen) 
end do 
Mag(k)=sqrt(Re*RetIm*Im) 
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print*,k,df*k,Mag(k) 
end do 
end 


4.5, MAXWELL’S EQUATIONS 


Thus far, we have dealt with the wave equation describing one 
scalar variable E, and thus we have dealt with an equation that 
involves second derivatives with respect to time. An alternative 
approach is to write two first-order equations. Although this can be 
done abstractly, for electromagnetic fields the curl equations of 
Maxwell introduce this approach with physical meaning for the 
variables involved. Since Yee’s pioneering paper [1], a number of 
applications of this method have been made. As a simple one- 
dimensional example, if the fields are uniform with respect to x and 
y, the curl equations reduce to equations involving first derivatives 
with respect to z and rf. One such field is described by 


z oy 4-3 

Seok 7 Bay (4-36) 
and 

on a 4-37 

ee 85 Be (4-37) 


Although with these continuous equations, we think of the two 
variables E, and H, as being defined at the same point in z and f, 
with the finite-difference equations this is not the case. In the first 
equation, for example, the derivative with respect to z is best 
thought of as being evaluated between the values of E,, and thus 
the values of H, are located in distance z between those values. 
The equation, however, locates the derivative of H, with respect to 
time 1 there. The derivative with respect to time of H,, should be 
evaluated in time between the values of H,, Similar considerations 
can be drawn from the finite- difference” version of the second 
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of of ° H 
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FIGURE 4-1. Mesh for Maxwell's equations in one spatial dimension. 


tion. The situation is summarized by Fig. 4-1. The finite-dif- 
nce approximations to (4-36) and (4-37) are 


E(k + 1,n) — E(k,n) 
h 


0.5, 5) — H(k + 05,2 - 05 
a eae + n+ 0.5) : ( n ) (4-38) 


H(k + 0.5,n + 0.5) — H(k — 0.5,n + 0.5) 
¥ h 
E(k,n + 1) — E(k,n) 


———— (4-39) 


first equation corresponds to the continuous equation at 
0.5,n)} and the second equation corresponds to (k,n + 0.5). 
al conditions can be imposed as values of E for n = 0 and as 
es of H for n = 0.5. 


The situation in two and three spatial equations follows the same 


ples, except that we must decide where to position the points 
hich the field variables are defined. As first shown by Yee, 
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FIGURE 4-2. Typical cell for TE wave 


these points can be so positioned that the finite-difference approxi- 
mations are central-difference approximations, which result in sec- 
ond-order accuracy. If a two-dimensional region is divided into 
square cells, Fig. 4-2 shows one arrangement of these points in a 
typical cell for the TE situation which involves E,, E,, and H., 
Although the relative position of the points at which the field 
variables are defined is fixed, their positions in the cell and thus 
with respect to the boundaries of the region can take several forms, 
For example, although the typical cell for the TM situation could 
take the form of a dual of Fig. 4-2, as shown in Fig. 4-3, a more 
convenient arrangement is shown in Fig. 4-4. The latter arrange- 
ment results in E_ being on the boundary, which is appropriate for 
a metallic-wall boundary. Just as in the one-spatial dimension case, 


FIGURE 4-3. One version of the typical cell for TM wave. 
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g A, e 
Hy Hy 
E. Ay E. 


FIGURE 4-4. Another version of the typical cell for TM wave. 


snetic field values are evaluated at times intermediate to the 
at which the electric field values are evaluated. 

he situation in three spatial equations can be arranged as 
yn in Fig. 4-5. The three electric field components and three 
netic field components are positioned to satisfy the require- 
ts for second-order approximations. In addition to their spatial 
ition, magnetic field values again are evaluated at times interme- 
ite to the times at which the electric field values are evaluated. 


IRE 4-5. Spatial arrangement of field components in three spatial 


nsions. 
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The finite-difference time-domain method requires a closed re- 
gion in order to have a finite number of nodes. Because problems 
involving open regions are important, much research has performed 
to derive special “absorbing boundary conditions” to be applied at 
the region boundaries to simulate open regions, This is still an 
active area of research. 


PROBLEM 


1. The finite-difference time-domain method is to be applied to the 
determination of f(x, 1), where 


rf 1 &F 
axe ar 
with the boundary conditions 
f(0)=0 and f(a)=0 


and the initial conditions 
f(x, 0) = sin rx/a. 
of 
a 9) =0. 


The spatial interval h is 


and the normalized parameter p is 
p= 9.5. 
In terms of the function F(k, ) defined by 
F(k,n) =f(kh,né), 
calculate F(2,2), 


REFERENCE 73 
R PROJECT 4-1 

purpose of this project is to examine the nature of waves in two 
nsions. Consider the TM electric field component E.(x, y,t), 


h satisfies the two-dimensional wave equation 


#E, 
ax? 


PE, 
ay 


1 @E, 
ce ar’ 


a square with sides L, with the initial conditions 


.t)=e"®, where R= ¥(x —L/4)° + (y - L/2y 
syst) 


dE 0)=0 
ay hs ) = 0. 


initial condition will propagate away from the vicinity of 
‘4, /2) and the resulting “pulse” observed at (L/2,L/2). 
an appropriate value for a, the pulse will be sharp enough 
Il enough localized) to be recognizable. Select (by experimenta- 
mn if necessary) a suitable value for a and evaluate the received 
ve E(L/2, L/2,1) for a long enough period of time to see both 
direct pulse and the first echo from the walls. Plot the response 
a function of the normalized time ct/L. What can you say about 
e shape of the received pulse? What about the echoes? 
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« STATIONARY CONDITIONS FOR FUNCTIONALS 


st common method for analyzing electromagnetic problems 
solving partial differential equations. Variational methods 
an alternative. A key concept in the variational calculus is that 
a functional. This concept is a generalization of that of a 
. Just as a function associates with each value of the 
endent variable(s) a numerical value, a functional associates 
| each function a numerical value. As a simple example, the 
of curve defined by the function f(x) can be expressed as 


d 2 
Lisa] = f° 1+(Z| de. (5-1) 
is is a special case of the functional 
d, 
If] = fens) ae. (5-2) 


= next key concept is that certain functions make the func- 
achieve a maximum or a minimum value. For example. the 
Potential # is known to minimize the electrostatic field energy. 
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The achievement of a maximum or minimum value, together with 
the equivalent of the saddle point (a function for which first-order 
variations of the functional are zero and yet the functional is 
neither maximum nor minimum), is referred to as stationarity. 

One method for determining the function that makes the func- 
tional stationary is to substitute for f the function f + ag, where g 
is a function that is arbitrary except that it satisfies homogeneous 
boundary conditions and @ is a parameter. Then f makes the 
functional stationary if and only if 


d 
qailf* ag]=0 fora=0. (5-3) 


As a simple one-dimensional example, to determine the minimum 
of 


dj 2 
aval -/(Z) de, with f(0)=0 and f(1)=1, (5-4) 


we replace f by f + ag, where g(0) = 0 and g(1) = 0. Then the 


integral can be written as 
is+ es] +1(f1 + 2¢{ (2) (3) ae] +71. (5-5) 


Use of (5-3) gives 


if df \(dg 
Malla)e-° a 
Integration by parts results in 
2 
Poe on 
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bracketed expression is zero because of the boundary condi- 
on g, and thus 


1 (af 
(Sr) eas =o. (5-8) 
this integral to be zero for arbitrary functions g, 
ar 
ae 3) 


more general functional in (5-2) can be analyzed in a similar 
ay and leads to the differential equation 
Fyf" + Fyf' + F, =0, (5-10) 


as the Euler-Lagrange equation. In a similar manner, many 
nals that are in the form of integrals are made stationary by 
ns that satisfy differential equations and this same approach 
be used to derive the differential equation. For example, the 
ional 
=" (fj +r dx, with f(0)=0 and f(1)=1 
alee ' 
(5-11) 


+ ag] =I[f] + zal f'(Z)(E) a+ fsea| +o7l[g]. 
(5-12) 


function f then can be shown to satisfy 


af 


oe (5-13) 
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5.2. EXAMPLE OF ELECTROSTATIC FIELD ENERGY 


An important example from electrostatics is the electrical field 
energy, 


mlocs rel = 5 ffel(se) + (2) + (2) ]ecera, 
(5-14) 


where the permittivity « = e(x,y,z) in general is a function of 
position. Use of the same method as in Section 5.1 results in 


W[d + ab] =1[4] +20f {fq | (2) p ($)(2] 
(3 


az lz | dedydz + a°I[W]. (5-15) 


Use of (5-3) again gives 


S14) ee) + (YF) + (2) aoe —o 


(5-16) 


This integral can be transformed by use of one of Green’s theorems 
(essentially integration by parts) to 


1S eel ax +3 : + Ee )]oaearae 0, 


(5-17) 


and thus 


[Stee] +2 [e2) +2 fett-0 
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in vector notation, 


V- (eV) = 0, (5-19) 


is of course the generalized form of Laplace’s equation that 


v-D=0. (5-20) 


. VARIATIONAL EXPRESSION FOR EIGENVALUES 


preceding examples of functionals involve a single integral. A 
nal for eigenvalues involves the ratio of two integrals. Con- 
first the one-dimensional eigenvalue equation 


+k*f=0 — with f(0)=0 and f(1)=0. (5-21) 
lication by f(x) and integration results in 
(ea Fac +h ['p ar = 0, (5-22) 
which 
fie aS te 
= % (5-23) 
ra 
0 
h integration by parts of the numerator, this can also be written 
d 2 
lz) « 
ee ‘ (5-24) 
[Pe 
0 
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Expressions (5-23) and (5-24) are functionals that depend on the 
function f(x) and are variational expressions. The functional ig 
minimized by the eigenfunction and the minimum value is the 
cigenvalue. Starting with (5-23) we consider as before 


1 af dg)? 
firs an (+0 dy 


K(f +ag) = — . (5-25) 


[+ aay a 
Tf we set 
a 2, 
dak UF + 28 )le-0 = 0, (5-26) 


the resulting equation can be simplified to 


ef a 
Gah Pae=s[ eee (5-27) 


which can be written in the form 


=0, (5-28) 


which with the notation of (5-23) is the same as (5-21). The 
minimization is of course only a local minimum and there are an 
infinite number of minima, one for each eigenvalue. 

In two dimensions, (5-23) and (5-24) generalize to expressions of 


the form 
[/ V7 fdxd) 
hn sf day (5-29) 


(Paces 
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i (Vf) - (VF) dedy ; (5-30) 


[f #? aeay 
ich correspond to the differential equation 
(V2 +k) f(x,y) =0. (5-31) 


In three dimensions, the generalizations are 


po — Lp fv fay de (5-32) 
[ff P aedyaz 


pe LLP ee aednde (5-33) 


i] ii fr dxdy dz 
hich correspond to the differential equation 


(V2 +k?) f(x, y,z) =0. (5-34) 


RAYLEIGH — RITZ METHOD 


Rayleigh—Ritz method is a general method to obtain an approx- 
to the function that makes a functional stationary. This 
jod is based on pioneering work by Lord Rayleigh in 1870 and 
ovements by Ritz in 1909. The function is approximated by a 
combination of known functions and the solution consists of 
mining the parameters in the combination. This determination 
substitution of the linear combination into the functional 

id then differentiation with respect to each parameter. 
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Conailbcraipins tie teatnpis e zero, as frequently is the case with eigenvalue problems, then 
js zero. The sum from (5-41) is substituted into the functional 
the derivative of the resulting expression with respect to each 
the parameters is set equal to zero, which yields a set of N 
ions. Solution of this set of equations yields the parameters 
thus the approximate solution for the function that makes the 


tional stationary. Substitution into (5-35) yields 


uni~ { (Ze) +8] 


Ee ff) ern 


it 


df 2 
w= f{(Z) “las with f(0)=0 and f(1) =1, 
(5-35) 


A function that satisfies the boundary conditions but does not 
minimize the functional is 


fol) =x. (5-36) 


A set of functions that can be used to improve the approximation is 
described by 


NON Af) 
fs) =x =x", 537) ry i Ean l'|(Z)( i) + Hednco] dv. (5-42) 


jrtkel dx }\ de 


These functions satisfy the homogeneous boundary values : : 
tiation with respect to cach a, gives 


Lal (Ze) + feyhteo] a 


=~ (BNE) +e 


f(0)=0 and f(1) =0, (5-38) 


and hence a linear combination of these functions can be added to 
the right side of (5-36) without changing the boundary values, The 


resulting approximation is dx, (5-43) 


F(#) =x + 0,(¥ — x?) +.a,(x — x3) + +++ ta(x— 9%"! 
i ) al ) ae is a set of linear algebraic equations from which the parame- 


oF G37) STS a, can be determined. 
N 
f(x) =x + Yaj(x-x'*), (5-40) ROBLEMS 
i=1 . 
A generalization of this expression is ‘Use the variational method on (5-1) to show that the shortest 
distance between two points in a plane is a straight line. 
N L * a 
S(%) = fo(x) + y a, f(x). (5-41) Use variational methods to show that the functional 
a=t 


where /,, satisfies the nonhomogeneous boundary conditions, the f,, 
for 1 <n <N satisfy homogeneous boundary conditions, and the 
a, are parameters to be determined. If the boundary conditions on 


df\2 
me tl fro 
with boundary conditions f(0)= 0 and /(1)=1, 
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tisfies the boundary conditions. Then evaluate the functional 
en by (5-24). Solve the differential equation (5-21) and deter- 

the true value of k7. What two conclusions would you draw 
om a comparison of these two values for k?? 


is made stationary by the solution of 


d? 
pad +f=0 (with the same boundary conditions). 


Note that this is not a minimum or maximum. Find the function 
f(x) by solution of the differential equation and evaluate the 


resulting value of the functional J. Evaluate the functional for 


the trial function R. Weinstock, Calculus of Variations (Reprint Edition), Dover Publica- 


New York, 1974. 
fila) =x 
and compare the result with the previous value. 
3. The functional J, defined by 


1 
If] = f (=) dx 
with the boundary conditions 
f(0)=0 and f(1)=1 
is made stationary by the function 
f(x) =x. 
Express as a polynomial in @ the functional for the function 
f(x) =x + a(x - x7). 


Sketch this polynomial. What conclusion you can draw from this 
sketch? 


4. For the eigenvalue equation (5-21) and the corresponding func- 
tional (5-24), consider the trial function 


f(x) =a + be + cx?, 


Choose the coefficients a, b, and ¢ so that the trial function 


IC CONCEPT OF FINITE ELEMENTS 


~element method developed initially in mechanical and 
ineering and most of the applications and, accordingly, 
d other publications have involved structures. The method 
ising importance in other fields, including electromagnet- 
al books should be of use to those interested in electro- 
tics [1-4]. 

basic concept of the finite-element method is that although 
or of a function may be complex when viewed over a 
gion, a simple approximation may suffice for a small subre- 
ie total region is divided into a number of nonoverlapping 
s called finite elements. In two dimensions we usually use 
and the simplest polygons are triangles and squares, 
6-1 shows a region divided into squares and Fig. 6-2 shows 
region divided into isoceles right triangles. Sometimes, as 
d in Fig. 6-3, a combination of triangles and squares is 
One of the advantages of using triangles is that a fairly 
region can be more easily approximately covered by a set 
les, as shown in Fig. 6-4. Regardless of the shape of the 
ts, the field is approximated by a different expression over 
lement, but where the edges of adjoining elements overlap, 
d representations must agree to maintain continuity of the 
- The equations to be solved are usually stated in terms not of 
Held variables but in terms of an integral-type functional such as 
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FIGURE 6-1. Division of a region into square elements. 


VAVAVAVA 
VAVAVAVA 
VAVAVAVA 
VAVAVAVA 


FIGURE 6-2. Division of a region into right isoceles triangles, 


FIGURE 6-3. Division of a region into triangles and squares. 
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FIGURE 6-4. Division of an irregular region into triangles. 


. The functional is chosen such that the field solution makes 
nctional stationary. The total functional is the sum of the 
al over each element. 


. FINITE ELEMENTS IN ONE DIMENSION 


ider the problem in one dimension of minimizing the func- 


Iif]= p\gfa (6-1) 


h the boundary conditions 


F(0) = 0 (6-2) 


fQ) = 1. (6-3) 


a manner similar to that of Chapter 2, divide the interval (0, 1) 
four subintervals and use the notation 


Fx = F(G.25k). (6-4) 
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If the function is approximated over the kth subinterval by 
F(X) = fan + ACh — fe-s)Le = 0.25(k - 1], (6-5) 


the derivative over this subinterval is approximately 


not very important, the functional values and the minimum value of 
he functional derived by the finite-element method happen to be 
act for this simple example. 

Now consider the slightly more complex example 


af 


2 
d ee ins | 2 
= mAh — Raid (6-6) i= {(z) “nila — 


Evaluation of the integrals over each of the subintervals yields us use the same approximation as before with four elements 
“ intervals). To the previously derived value for the first func- 
If] =4(f, ‘. +4f2-fi) +4 -he y +41 fry. (6-7) we must add the result of integrating the second term, which 
an be evaluated as 


To find the functional values that minimize this expression, differ- 
entiation with respect to each of these values gives 


[ra = n(2fF + Sif + 262 + fofs + 265 + fa +1). (6-15) 


dl 
a Bf) + Bf — fa), (6-8) total functional then is 
dl 
de Ah — fi) + 82 fa), (6) I[f] = 8.17f? — 7.92, fa + 8.172 - 7.92f af 
and + 8.17f2 — 7.925 + 4.08. (6-16) 
dl 
a 8(fs — fa) + 85 — 1). (6-10) ting the derivatives of the functional with respect to f,, fy, and 


qual to zero leads to the equations 


16.34f, — 7.92f, = 0 


f, = 0.25, (6-11) ~7.92f, + 16.34f, — 7.92f, = 0 
z= 0. 12) 

f= 05, (oa) —7.92f + 16.34f; = 7.92, (6-17) 
and 

f; = 0.75. (6-13) hich have the solutions 
Substitution of these values into (6-7) gives the value of 1 for /[f} fy =0.21 
Two observations can be made about this solution. First the equa- f,=0.44 
tions derived by minimizing the approximate expression for the lates 
functional can be seen to be the same equations that result from = 06, (6-18) 


application of the finite-difference method to the differential equa- 
tion whose solution minimizes the functional, Second and probably 


vhich agree well with the analytic solution. 
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6.3. LINEAR INTERPOLATION FOR ISOCELES 
RIGHT TRIANGLES 


In two dimensions the simplest element is the isoceles right trian- 
gle. Consider the typical such triangle shown in Fig. 6-5, For 
convenience while analyzing this triangle, the origin of the coordj- 
nates is taken at one vertex. The simplest interpolation to use is the 
first-order formula 

o(x,y)=at+be+cy. (6-19) 
There are three unknown parameters (a, b, and c), so it is natural 
to specify the interpolation in terms of the value of the function at 
each of the three vertices of the triangle. Solution for the parame- 
ters results in 


o,— % 
Rta 


(x, 9) = by + “y, (6-20) 


which can be written in the form 


Hod (0-§ Four fove(E}on can 


usually referred to as a nodal expansion. When the field quantity 
being sought is the electrostatic potential, the electric field energy is 
a natural functional to use because the potential is known to 
minimize this energy. This energy in two dimensions for uniform 


O 1 * 
FIGURE 6-5. A typical right isoceles triangle. 
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5 |[ (Vs) - (Vo) dedy, (6-22) 


(6-23) 


-$|(2) + (Fl Jao 


energy for the isoceles right triangle is found by substitution of 
(6-20) into (6.23), which yields 


ba 2 F3 

W = F[(4o - 41)" + (4p ~ 42)']- (6-24) 
expression clearly does not depend on the orientation of the 
le and applies to other orientations. If the total region is 
led into a number of right isoceles triangles, the energy for 
element is proportional to the sum of squares of the differ- 
s of the vertex values, using the general form of (6-24), The 
energy then is found by addition of the energy contributions 
mm each element. 

The true potential minimizes the total energy, and therefore the 
pproximate solution is determined by differentiation with respect 
© the parameters, which in this case are the unconstrained values 


AVA 
WA 


FIGURE 6-6. A typical interior node with triangular elements. 
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of the potential at the triangle vertices. Note that the boundary print* 
conditions constrain values at the nodes that fall on the boundary, print*,' what is a/d?' 
The energy is the sum of a number of terms, some of which involve read*, add 
@(i, j) and some of which do not. For a typical interior point (x, y) print* 
as illustrated in Fig. 6-6, print*,'what is b/d?" 
read*,bdd 
‘ 8 a 2 rint* 
oneray terms that involved — leo 7 41] zl - $2] Bet atk Shou many values of M2?" 
€ a. © . read*,WW 
+5140 - 43] + 5140-4] . print® 
print *,' enter tolerance’ 
(6-25) read*,tol 
print* 
Setting the derivative with respect to &, equal to zero yields print * 
R=1.4 
_ bi + br + b3 + bs print parameters 
by = 4 : (6-26) print*,' shielded capacitor’ 
print*,' a/d is',add 
which is the same equation derived earlier with the finite-difference print*,' b/d is',bdd 
method. The following FORTRAN program uses the finite-element print*,' tolerance is',tol 
method with isoceles right triangles to determine the shielded print*,'' 
parallel-plate transmission line capacitance problem considered in print*,' M original extrapolations' 
Chapter 2. print*,' 
start calculation 
program shcapfe M=1 
C capacitance of shielded capacitor W=1 
c calculated with finite elements do while (W.LE.WW) 
c Line integral of E field nx=add*M 
ie solution by iteration ny=bdd*M 
real add,bdd,R,delta,cap,old,z set initial and boundary conditions 
integer C,W,WW,M,nx,ny,i,j,k,kk do j=0,ny 
real V(200,200),ZZ(10,10) do i=0,nx 
c input parameters VCi,j)=0 
print* end do 
print*® end do 
printx,' shielded capacitor’ do j=1,ny-1 
print*,' finite-element approximation’ do i1=0,nx-1 
prints,’ capacitance from Line integral’ VCi,j)=1 
prints,' solution by iteration’ end do 
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end do 
cap=1 
old=0 
k=1 
kk=1 
do while (abs(old-cap).GT.tol) 
old=cap 
do while (k.LE.kk) 
do j=1,ny-1 
do i=0,nx-1 
if (j.EQ.M) then 
if (i.GT.M) then 
2=VC4-1,j) 4VCi41,)) 
z=z+VCi,j-1)4V C4, 541) 
delta=z/4-V(i,j) 
VCi,j)=VCi,j)+R*delta 
end if 
else 
if (i.E@.0) then 
zZ=2*V(14+1,j)+VCi,j-1) 
z=24VCi,j+1) 
delta=z/4-V(i,j) 
VCi,j)=VCi,j)+Redelta 
else 
2=VC4-1,j)4VCi41,5) 
Z=z24VCi,j-1)4+VCi,j+1) 
delta=z/4-V(i,j) 
VCi,j=VCi,j)+Redelta 
end if 
end if 
end do 
end do 
k=k+1 
end do 
Cc start calculation of capacitance 
z=0 
do j=0,ny-1 
do i=0,nx-1 
z=2+(VCi41,j)-VGi,j))* 
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(Git1,j)-VGi,j)) 
zezt(VCit1, j+1)-VCG,jt1))* 
(VCi41, 541) -VGG, 5410) 
z=2t(VCi,jt1)-VCI,j))* 
(VCi,jt1)-VCG,j)) 
z=zt(VCi+1,jt1)-VCit1,j))* 
(Vv Ci41,5+1)-VCI41,5)) 
end do 
end do 
cap=8.854187*z/2. 
kk=2*kk 
end do 
ZZ(W,1)=cap 
extrapolation calculation 
c=1 
do j=2,W 
C=2*C 
ZZ(W, j)=(C#ZZ CW, j-1 -ZZ0W-1, 5-1 / 00-1) 
end do 
print results 
print*,M,ZZ(W,1),ZZ(W,W) 
W=W+1 
M=2*™ 
end do 
end 


An interesting situation arises at boundary points for which the 
undary condition is that the normal derivative is zero. For a point 
unknown value on a boundary as illustrated in Fig. 6-7, 
energy terms that involve 
€ 2 € 2 € 2 
= Gb — by) + Fb ~ ba) + 5 (ba be). (6-27) 
tting the derivative with respect to (x, y) equal to zero yields 


_ dit dbs + 265 


bo 4 


(6-28) 
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3 0 1 Boundary ¥ 
3) 2 
4 
FIGURE 6-7. A node with unknown value on a boundary. 
te) pi x 


Note that although the boundary condition that the normal deriva- 
tive is zero was not used, the result is the same as the finite-dif- 
ference result, where we used that boundary condition explicitly, 


FIGURE 6-9. A square element divided into triangles. 


erpolation over the square. Such an interpolation that uses the 
r nodal values should have four coefficients. One such approach 
6.4. SQUARE ELEMENTS bilinear interpolation 
Consider now the square element shown in Fig. 6-8. One simple H(x,y) =a + bx + cy + dey. (6-30) 
approach to determining an approximation to the energy of this 
element is to divide the square into two triangles, as indicated in 
Fig. 6-9. Addition of the energy for the two subelements gives 


Evaluation of the four coefficients and rearranging as a nodal 
yansion gives 


ste =(1-F)lt— a) Calta) 


Clabes AlR}o> om 


ation of the energy using this expression gives 


= S 2 2 2 2} 
W = = [(by ~ 41)? + (bs ~ bo)? + (bq ~ 43) + (4 ~ 42)']- 
(6-29) 


This approach, although workable, does not correspond to a smooth 


y ag 2 2 
4 = el(40- 41) + (bo — b:)(b3 — b2) + ($3 — 2) 
: . + (by — bs)" + (by ~ $3), ~ 42) + (6, ~ 2) | (6-32) 
derivative of W with respect to one nodal value #y can be 
ated as 
i) 1 x a. E 
sa = Zl4b0 — b) — 262 — 43). (6-33) 
FIGURE 6-8. A typical square element. 9b 6 
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4 3 2 Oe.) 
(2. ¥2) @) 
@® (*3.¥3) 
6 7 8 FIGURE 6-11. The general triangular clement. 


FIGURE 6-10. A typical interior node with square elements. 


esentation to agree with the functional values at the vertices 


Using a similar expression for the other three squares shown in 


Fig. 6-10, the derivative of the energy for the four squares is 
at+bx, +o, =4, 


awe at bx, + cy, = 
Ibe = G (loo — 26, — 2) — 2, — 2d, 
a a + bx; + cy; = 43. (6-36) 
—2bs — 2b4 — 2h, — 24s). (6-34) & a a 
Nn matrix notation, 
Since ¢ is selected to minimize the total energy, the left side of 
; fon i 1 x, y,//a@ co 
this expression is zero and 
1X2 y2))b) =] $2). (6-37) 
1 x; y3]}e¢ bs 


4 _ bi + bo + b3 + by + bs + be + b7 + Os 
a= . 


(6-35) 


e e determinant of the coefficients is 
lox yy 
6.5. GENERAL TRIANGULAR ELEMENTS 1 x, ¥2|/=2A, (6-38) 
1 ox; ¥3 


Next we generalize the situation from the isoceles right triangle t0 
the general triangle, as illustrated in Fig. 6-11. An excellent account 
of the required analysis is given in Reference 1. We would like to 
use the same linear representation (6-19) that was used with the 
isoceles right triangle. The equations that must be satisfied for this 


e A is the area of the triangle. Except for the special case in 
h the three vertices are collinear and the area is zero, the 
uations can be solved for a, b, and c. Then (x, y) can be 
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expressed in matrix form as 


o, 
(x.y) =[%, w2 Hs] $2], 
$3 
where the w-matrix is given by 
L Ry ik 
[¥) w ws) =[1 x y]]l om 
1 x3; 3 


The expression for (x, y) is equivalent to 
3 
O(x,¥) = L(x, v)be- 
k=1 


Evaluation of the yy, gives 
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nto (6-22) to yield 
€ 
W=~L 5b Sy, 6-45) 
(6-39) 2 x Pe jh ( 
in = ff (Vu,) = (Wu,) deay. (6-46) 


matrix coefficients S,, depend only on the shape of the 
ngle, not on its size or location and orientation with respect to 
coordinate system. For example, 
2= ff aby Oh 4, A 

“Ox ax 


6-47 
sae ray" (6-47) 


| dx dy, 
(6-41) 


is 


1 
Sp= Fallin —Ys)(¥s — Yi) + (%3 — ¥2)(41 —x;)], (6-48) 


1 
UCase ¥) = 55 [(¥a¥s — ¥a¥a) + (va — ya) + (3 — 42) 9] 
ich can be evaluated (see Reference 1) as 
oly) = soll — V3) + (ys — ya + (4, —45)¥ 
2! | eral 1 3) ( 3 ) ( 1 3) ] Sum —tcotd,, (6-49) 
U(x, ¥) = sls - ¥294) + (Y —y2)e + (42 — x) ¥I- ere 0; is the interior angle at vertex 3 of the triangle. The other 


The ws, are interpolation functions in the sense that 
é 1 iff=k 

Xr Vy) = 

eM Ny tga, 


To evaluate the energy integral, substitute 


Vb(x.y) = LON, 
k 


(6-42) coefficients can be evaluated in a similar way, which leads to 
1 cot 0, + cot 0; —cot 0; —cot 6, 
ae —cot 6; cot @, + cot 4; —cot 6, 
—cot 6, —cot 4, cot @, + cot 4, 


(ng (6-50) 


The energy calculated for one triangle from (6-45) is a quadratic 
ion of the potential values at the vertices. The total energy for 
entire region is the sum of the energies of each triangular 
ment and hence is a quadratic function of all the potential 


(6-44) 
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values. The minimum energy is found by differentiation of this 
function, which results in a linear combination of potential values 
for each node. Thus the result is a set of linear equations that can 
be solved for the unknown potential values. Because the matrix of 
this set of equations is sparse, iterative methods of solutions nor- 
mally are used, 


6.6. HIGHER-ORDER INTERPOLATION WITH TRIANGLES 


We wish to consider use of higher-order polynomials, but with our 
current approach it is clear that the algebraic computations will be 
complex. Clearly, we need a more systematic way to write the 
equations. The method we will follow involves the use of simplex 
coordinates, sometimes referred to as area coordinates. In two 
dimensions we normally define a point by the two spatial coordi- 
nates x and y. With simplex coordinates, we specify a point by 
three coordinates s,, s,, and s;. These coordinates are defined 
locally with respect to one triangle and we use a different set to 
define a point in another triangle. Figure 6-12 shows a typical 
triangle and the s; coordinates. Each of the three coordinates is 
defined as the ratio of the distance of the point from a triangle side 
to the distance of the opposite vertex from that side. Defined in this 
way it is clear that each s, varies from 0 to 1 as the point moves 
from the side opposite vertex j to the vertex j. The redundancy 


3 
FIGURE 6-12. Simplex coordinates for the triangle. 
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5, +5,+5,=1. (6-51) 


cause each simplex coordinate is the ratio of the altitudes of two 
ngles with the same base, it clearly is the ratio of the areas (from 
vyhich comes the other name, area coordinate). For example, s, can 
expressed as 


1 #& ¥ 
1 ox. 2 
1 x3 Ys 
= 6-52 
ay I (6-52) 
Tox, ¥2 
1 x; Ys 


Rather than use this expression to evaluate the simplex coordi- 
Mate s, as a function of x and y, we can see from Fig. 6-12 that if 
point P is at vertex 2 or 3, s, is zero and if the point P is at 
ertex 1, s, is unity. Similar considerations apply to the other two 
nplex coordinates, s, and s,. Therefore, the simplex coordinates 
equal to the interpolation functions defined by (6-42). From 


wy = 5, 
Wy = $y 
W; = S3. (6-53) 


is suggests that higher-order interpolation formulas could be 
lefined in terms of the simplex coordinates. Consider the quadratic 
‘Mterpolation expression 


(x,y) =a + bx + cy + dx? + exy + fy’. (6-54) 
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3 
FIGURE 6-13. Nodes for the quadratic representation 


Because there are six coefficients, we need six nodes, and Fig. 6-13 
shows a logical way to locate these nodes. Interpolation functions 
for these nodes can be constructed from expressions that are 
quadratic in the simplex coordinates. The quadratic interpolation 
function #,, for example, must be unity for s, equal to unity and 
zero for s, equal to zero and one-half. Hence 

Y, =5\(2s, — 1). (6-55) 
The other five interpolation functions can be constructed in the 
same way and the result is 


Wy, = 5,(28, — 1) 
Wz = $2(25, — 1) 
% = 53(253 — 1) 


(6-56) 
hy = 45,82 
Ws = 45255 
Ue = 45,53. 


The next higher-degree polynomial is the cubic, which has 10 terms. 
There are not 10 logical nodes on the triangle edges, so a node 
inside the triangle is used, as shown in Fig. 6-14. Similar functions 
can be constructed from the simplex coordinates for cubic and 
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9.7. NODAL EXPANSIONS AND THE WEAK FORMULATION 


have written several local approximations in the form of nodal 
sions. We can write one global expansion in the form 


N 
(x,y) = Wooley) + Z etal, 9), (6-57) 


€ w(x, y) satisfies the boundary conditions, N is the total 
lumber of unknown nodal values, @, are these nodal values, and 
X, y) are the interpolation formulas, which are derived from the 
formulas. This global expansion can be used with the 
igh—Ritz approach to select the #, that minimize the total 
gy. This is simply another way of describing the approach that 
have followed before. Use of the nodal expansion as an approxi- 
ion to the weak formulation. 


J Vo - Vis dx dy (6-58) 
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results in 


N 
Zoe ff Vs “Wy, dedy = — [[ Vig Vu, dedy. (6.59) 


6.8. TIME-DEPENDENT VARIABLES 
The finite-element situations considered up to now have been 
static. One way to handle dynamic situations is to use finite-element 
expressions for spatial variation but with nodal variables that de- 
pend on time. For the wave equation 

PE 4 ; 

ae 7 ¢v’E with E=Oonthe boundary. (6-60) 
Substitution of an approximation by the nodal expansion 


Blt, ¥ 201) = LER Wale ¥>2) (6-61) 
gives 
LEI, ae CLE. (6-62) 
Multiplication by #, and integration results in 
LE; J ff vit de dy dz = PLE sf [[ON7e. (6-63) 
which can be written through one of Green’s theorems, 
DEL fff tle dedy de = ~c? DE, fff Vi, Ve dedyde. (664) 


This expression is a weak formulation and has the advantage that 
only first spatial derivatives of the expansion functions are involved. 
Using a finite-difference approach with the time variable is one of 
several ways to solve this set of ordinary differential equations. 
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LEMS 


the element that is the nonisoceles right triangle shown 
, calculate the energy in terms of the nodal values. 


te 


0 hy oy 
PROBLEM 6-1. Right triangular element. 


For the square element shown with the nodal potential values in 
the following figure, calculate the energy for the two methods 
discussed in the text. 


o=3 o=6 
PROBLEM 6-2. Square element. 


Calculate the energy of the nonsquare rectangular element shown 
below. 


3 2 
hy 

6 
#0 hy i 


PROBLEM 6-3. Rectangular element. 
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4. The finite-element method is to be used with square elements 
and bilinear interpolation to calculate electrostatic potential, 
Derive the equation to be used to determine a nodal value when 
the node is on the boundary between two media with different 
permittivities. 


P. Silvester and R. L. Ferrari, Finite Elements for Electrical Engi- 
, Cambridge University Press, New York, 1983. 

_ C. W. Steele, Numerical Computation of Electric and Magnetic Fields, 
Van Nostrand Reinhold, New York, 1987, 

_ A. J. Davies, The Finite Element Method: A First Approach, Oxford 

University Press, Oxford, 1980. 

B. Bickford, A First Course in the Finite Element Method, Richard 

COMPUTER PROJECT 6-1 Irwin, Homewood, IIl., 1990. 

For the shielded microstrip shown in the figure, use the finite-ele- 

ment method to calculate the capacitance per unit length. Use 

square elements with bilinear interpolation functions. Express the 

result in pF/m to the nearest 0.1 pF/m. Take advantage of the 

symmetry. The assumptions are: 


1. €, = 95 €g. 


2. The microstrip (conducting strip on the dielectric material) is 
very thin. 


3. The conductors are lossless. 
4. The dielectric material is linear, lossless, and isotropic. 


0.4.cm 


COMPUTER PROJECT 6-1. Cross section of shielded microstrip transmis- 
sion line. 


1. LINEAR OPERATORS 
tromagnetic problems usually involve solution of linear partial 


ntial equations or integral equations. The general form of 
a linear equation is the operator equation 


f(R) = 8(R), (7-1) 


is a linear operator, f is the function to be calculated, 
g is a known function. Two electrostatic examples are 


otk) 


0 


V°4(R) = — (7-2) 


the operator involves derivatives and is referred to as a 
tial operator, and 


DI seca) PR = OR) (73) 


© the operator is called an integral operator. The general 
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integral operator can be written as DETERMINATION OF THE PARAMETERS 

e two sides of (7-6) cannot in general be equal everywhere. 
parameters @, are to be determined so that the resulting linear 
abination is a good approximation in some sense to the desired 
munction f(x). We can impose at most N constraints from which 
a, can be calculated. One approach, called point matching, sets 
two sides of (7-6) equal at N points R,, giving the set of 


ations 


S{{CCR RFR) BR’ = 8(R), (7-4) 


where the kernel G(R,R’) is known as a Green's function. 


7.2, APPROXIMATION BY EXPANSION IN BASIS FUNCTIONS 


To solve problems such as those described by (7-1), the method of N 

moments [1] begins by approximating the unknown function by a X a,g,(Rj) =A(R,)- (7-9) 
linear combination of known functions u(R) in the form k=3 

is a set of linear algebraic equations that can be solved to 
obtain the parameters @,. Substitution of the resulting set of 
meters into (7-5) completes determination of the approximate 
where the functions Y are known as basis or expansion functions olution. 

and are so selected that with appropriate values for the parameters _A more general approach is to multiply each side of (7-6) by one 
a,, the right side of this equation is a reasonably accurate approxi- of a set of weighting functions w, and integrate over the domain of 
mation to the left side. The interpolation functions in Chapter 6 are R to yield 

special cases of basis functions, The function w, satisfies the same 
boundary conditions as the function f but otherwise can be se- 
lected arbitrarily, and the basis functions satisfy homogeneous 
boundary conditions. Substitution of (7-5) into (7-1) and use of the 
linearity of the operator results in an equation of the form 


N 
FR) = ho(R) + x a(R), (75) 
tl 


13 ay ff fwj(R)sa(R) R= ff fw,(R)h(R) dR (7-10) 
=1 


the set of algebraic equations to be solved for the parameters. 


N point-matching method can be considered to be a special case 
¥ a,g,(R) = A(R), (7-6) this method by taking the weighting functions to be 
k=l 
hie w)(R) = 6(R — R,). (7-11) 
g,(R) =w,(R) (7-7) important approach, known as Galerkin’s method, uses as 
ata hting functions the same set of functions as used for the basis 
A(R) = g(R) — -7(R).- (78) _ 


There are two fundamental classes of basis functions, entire-domain w)(R) = ¥)(R). (7-12) 
functions, which span the entire domain, and subdomain functions, 


each of which is zero everywhere except in a subdomain. N mathematical circles, the name method of weighted residuals is 
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frequently used, because if the residual Res(R) is defined as 
N 
Res(R) = xz @,8,(R) — h(R), (7-13) 
=1 
the integrated residual weighted with w, is set equal to zero, so that 


f f f w)(R)Res(R) dR = 0, (7-14) 


which clearly is equivalent to (7-10). The set of equations (7-10) can 
be written in matrix form as 


Aa =B, (7-15) 


where 


Ay = [ff (R)g.(R) BR (7-16) 


and 


B,= ff fw(R)A(R) &R. (7-17) 


7.4, DIFFERENTIAL OPERATORS 


Solution of problems involving differential operators corresponds to 
the solution of differential equations. As a simple one-dimensional 
example, consider 


2 
ql) = #4, (7-18) 
with the boundary conditions 
f(0)=0 and f(1)=0. (7-19) 


These boundary conditions results in y/, being zero. Basis functions 
with continuous derivatives at least to the first order are desirable 
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avoid delta functions. A suitable family of basis functions for this 
roblem is 


u(x) =x—x**) fork =1,2,...,N (7-20) 


1d using the Galerkin approach, we take the weighting functions 
the same functions 


w(x) =x —xI7h, (7-21) 
, for example, N = 2 gives the approximation 
f(x) = a(x — x?) + a(x — x3). (7-22) 
itution into (7-9) yields for the matrix coefficients 
Ay = fe -)Q)de= 5 


Ayy= ['(x—x°)(6x) de = 3 


nie 


Ay = ['(x-2°)Q2) dr = 
° (7-23) 

Az = '(x—x)(6x) de = 43 

B= [xe -2)8) de = 5 


B, = f'(x— 2x) de= 


fhe matrix equation for the parameters thus is 


Wie tal 


| Ie] 7 [: } (7-24) 


118 = METHOD OF MOMENTS 


whose solution yields 


—— 

al | 
sl- 

eas 


(7-25) 


and thus the approximation is 


F(x) = —yo(4 — x?) + g(x — x4). (7-26) 


7.5. INTEGRAL OPERATORS 


The second class of operators, the integral operators, involve 
weighted integrals of the unknown function in the form 


a(R) = [[[G(R, RVR) &R, (7-27) 


For this operator (7-8) becomes 


8(R) = [[fG(R RY (RY) aR (7-28) 


The resulting matrix coefficients are 


An= {fff [fu GRR (R) GR ER 


and 
B= J i f w;(R)h(R) dR. 
With the Galerkin method these matrix coefficients are 


An=ffff [fy RGRRW(R)ER AR — (731) 


and 


B, = ff [uj(R)n(R) aR. (7-32) 
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7.6. PULSE FUNCTIONS 


of the most widely used type of subdomain function is the unit 
se function, which is equal to unity over one subdomain and zero 
here. If the typical subdomain is denoted by S,, 


1 ifRisin S, 

R)= 7-33 

He(R) (5 elsewhere ( ) 

Sulse functions are especially useful with integral operators, where 

they frequently are used for both basis and weighting functions. 
With pulse basis functions, (7-28) becomes 


8(R) = fff GRR) aR, (7-34) 


and with pulse functions for both basis and weighting functions, the 
trix coefficients given by (7-31) and (7-32) become 


An= ff fff [ORR aR w&R 


(7-35) 


(7-36) 


B, = ff, rR) arr. 


ause with this approach we have divided the source region into 
t of “finite elements” with simple approximations to the un- 

n function over each element, this approach sometimes is 
sidered to be a part of the finite-element method. The differ- 
between the present approach and the finite element method 
ssed earlier is that here we usually are dividing just the source 
region and not the entire domain. 


PARALLEL-PLATE CAPACITOR IN TWO DIMENSIONS 


‘© now consider a two-dimensional capacitor that is the cross 
Section of a parallel-plate transmission line as shown by Fig. 7-1. 
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ae 


FIGURE 7-1. Cross section of paralle! plate transmission line. 


Denote the width by W and the separation by D. In two dimen- 
sions, the relation equivalent to (7-3) is 


1 
(x,y) = 5p J fol’ yin (x—-x) +(y—y) ae’. 
0 
(737) 


Take the potential on the plate at y’=D/2 as +1 and the 
potential on the plate at y’ = —D/2 as —1. We wish to determine 
the capacitance per unit length of this capacitor by first solving the 
integral equation (7-37) to determine the charge densities on the 
plates and then integrating the densities to determine the total 
charges, from which the capacitance can be determined. If the 
plates are assumed to be very thin, the integration is over the 
one-dimensional region of the plates and can be written as two line 
integrals, one for y’ = D/2 and one for y’ = —D/2, 


[Om plas D/2)in Vix —xY + (Y — D/2Y ee” 


~0.51 


és.) = = 


27rEy 
i 


2arey /-osw 


(7-38) 


[" pat, -D/2)in Ve — HY (y + D/P 
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where now the charge density p, corresponds to a surface charge 
density. Because 
p(x’, -D/2) = —p,(x', D/2), (739) 


(7-38) can be written as 


1 posw (x =x")? + (y + D/2P° 
s(x", D/2)h i) H, 
se vay (ea Py De 


27e, 
(7-40) 
ause of the symmetry with respect to x, 
p,(—x", D/2) = p,(x", D/2) (7-41) 
equation can be written as 
osWw 
(x.y) =f p(x’, D/2)G(x, ¥.x')dx’. (7-42) 


, 2 2 
(2, 95<)= 1 o (x am a a) 
27 \ @oay)+(y - bpp 
, 2 2 
an (x' +x) + (y + D/2) _ (7-43) 
) (+x) + (yD 


ause the potential on the upper plate is unity, we can set 
= D/2 and get 


OSW 
1=f ON p(x )G(x, x) de’ for0 <x <0.5W (7-44) 
0 
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and G can be written 


i (x! —x)}* + D? (x +x) + D? 
G(x,x') = Tink; / @—xp + of (x +a) L 


(7-45) 


where we have simplified the notation for p, and G. Equations 
(7-44) and (7-45) describe the specific equation to be solved. 

One approach to the solution of the integral equation is to use 
pulse functions for the basis functions and point matching at the 
midpoints of the subintervals to derive the equations for the param- 
eters. First, divide the interval (0,W/2) into N subintervals of 
length 


h is 7-46 
he a (7-46) 
and then the kth subinterval is 

S, = (kh —h, kh). (7-47) 


An approximation in pulse functions leads to 
N 
p(x) = La, P,(x), (7-48) 
k=l 


where 


1 ifkh—-h<x<kh 


(7-49) 
elsewhere. 


x) ={ 


If (7-48) is substituted into (7-42) and if the two sides are point 
matched at the midpoint, 


y= ih — 0.5h, (7-50) 
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the equation becomes 


N, (751) 


yr kh 
Laf Ga,,.x)de'=1 for j = 1,2, 
k=1 0 Tkh-h 
yhere from (7-45), 


1 (x! —x,) + D? (x +x) + D? 
x;,x') = In 7 — +p = 
27€, ('—% (xe +x,)° 


(7-52) 


relations (7-51) and (7-52) define the set of algebraic equations 
must be solved to determine the a, parameters. Because the 
ntial difference between the two plates is 2, the capacitance, 
h is 


C= i (7-53) 
becomes 
c=8f"" plxyar’ (7-54) 
—w/2 
and because of the symmetry 
w/2 
C= fp (x) de’, Wa) 
o 


N 
C= Vagh. (7-56) 
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It is convenient to rewrite (7-51) as 


yr haf G(x, x)ae =4 
a,h— - = 
fe” hr Hae 


and then take a,h as the unknown parameters. Because a, is the 
charge density in the kth subdomain, a@,/ is approximately the 
charge in that subdomain. Then if (7-57) is written in matrix form as 


Aah = 1, 
the matrix coefficients are 
1 kn 
Ay ah, Ge ) de 
and 
B,=1. 


Substitution of (7-52) into (7-59) and use of (7-50) leads to 
(x' — jh + 0.5K)? + D? 


kh 
™ 3 
khoh (x' — jh + 0.5h) 


7 x'+jh —0.5h) + D? 
fin yi es eat — de 
Kho (x' + jh — 0.Sh) 


With the variable change 


1 1 
- 2re,|h 


Ay 


1 
h 


‘ 


x’ —jh +0.5h 
D 


in the first integral, and in the second integral 


x + jh —0.5h 
a= D 


> 
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(7-61) becomes 
(7-57) 1 k-7+0.5) pee 
Ay = —— In s— a; 
ik Qmegy ft -j-05) A B 
< 
yk+j-05), | B+ P 
— In dB|, 7-64 
| ae B? B eet 
(7-58) 
here 
h W/D 
(7-59) = ao (7-65) 


(7-60) ition (7-64) can be expressed as 


2reEjy 


i [Filyk — yi + O.5y) — F\(yk — yi — 05y) 
+F\(yk + yj — 0.Sy) — Fy(yk + yi — 15y)]. (7-66) 


ere F, is an indefinite integral 


. (7-61) B+P 
F(B) = [my —ar— ab, (7-67) 
which can be evaluated as 
(7-62) Pee 
F(B) = Bin —— + tan~! B. (7-68) 


‘or an alternative solution, the Galerkin method uses for the 
testing function w, the same function as the basis function. With 
Pulse functions for both basis and testing function, instead of (7-51) 


(7-63) 
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the basic equation becomes 


N ‘t 
Ya,f” 


fr G(x, x!) de’ de = J” dx for j= 1,2 
k=) “Uhh? kh—-h hah »2,-..,N, 


(7-69) 


where G is given by (7-45). The right side of the equation is equal 
to A and thus the equation can be written as 


jh 
Las ah (Glee) de de= 1 forj=1,2,...,N, 
=I 
(7-70) 


If as before the unknowns are taken as @,h, the matrix coefficients 
are 


1 
An= 7 in Se {Gla x?) de’ de (7-71) 
and 


B= 1. (7-72) 


By making changes of variable similar to (7-62) and (7-63), the right 
side of (7-71) can be evaluated as 


Daregy? Falk ~ yj + ¥) — 2Fi(yk — yi) 


Ay = 
+F,(yk — yj — y) + F,(yk + yi) 


—2F,(yk + yj — y) + F,(yk + yi — 2y)], (7-73) 
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ere F,, which is an indefinite integral of F,, is given by 
0 if B =0 
| eB e . 
—— VBP + +1 ~ > Inyp? + Btan-'p if p #0. 
(7-74) 


following table shows the results of using these two methods to 
te the capacitance for the particular case W = D. 


Capacitance (pF /m) 
N Point Matching Pulse Weighting 
1 17.1098 17.7084 
2 17.8677 18.1861 
4 18.2858 18.4444 
8 18.5059 18.5851 
16 18.6188 18.6583 


I e correct value can be shown by further computation and extrap- 
tion to be approximately 18.73 pF/m. Figure 7-2 shows the 


rge distribution for pulse weighting with N = 16. 


|—Graph A 


0 0.25 05 0.75 af 
x/0.5W 


FIGURE 7-2. Charge distribution. 
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PROBLEMS 


1. Use the method of moments to solve 


da 
mls) = -—6x, 
with the boundary conditions 
f(0)=0 and f(1)=0. 
Use N = 1 and take the one basis function as 
f(x) =x —x?. 
Evaluate a, for the two weighting Pieter 


w, =x —x? 
and 
w,=1. 


By integration, the true solution is 


f(x) =x-x. 


Plot the true solution and the two approximate solutions. 


2. Use the method of moments to solve 


2 
gills) = 27, 


with the boundary conditions 


f0)=0 and f(1)=0. 


Use the Galerkin method with the two basis functions 


f(a) =* = 3? 
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and 
f(x) =x — 23. 


Compare the values for x = 0.5 of the true solution and the 
approximate solution, 


COMPUTER PROJECT 7-1 


‘Calculate the capacitance per unit length of the unshielded paral- 
Jel-plate transmission line using the Galerkin method with pulse 
tions and the parameter ratio 


one basis function 


w(x) = 


For the second solution, evaluate the required integrals numeri- 
cally. Compare the two solution methods. 


i 
COMPUTER PROJECT 7-2 


culate the capacitance with respect to infinity of a thin flat 
Square plate with 1-m sides. Divide the square into M* squares and 
‘use the moment method with two-dimensional pulse functions for 
basis functions. Do point matching at the center of cach small 
‘Square. Evaluate the integrals involving the Green's function nu- 
_Merically. 


_ REFERENCE 


1. R. F. Harrington, Field Computation by Moment Methods (Reprint 
Edition), R. E. Krieger, Malabar, Fla., 1987; Original Edition, 1968. 


1. FUNDAMENTAL SCATTERING EQUATIONS 


he solution of scattering problems is one of the major applications 
the method of moments. In such a problem, a wave called the 
ident wave strikes an object, which we will usually assume is a 
perfectly conducting body. This wave causes a current to flow on 
he object’s surface, and in turn, this current causes a wave called 
scattered wave, which can be expressed as 


E(R) = [[G(R.R)I(R) aR’, (8-1) 


here the integration is over the surface of the body and J, is the 
ace current density (in A/m). If the body is a perfect conductor, 
e tangential value of the total E field is zero on the surface so 


Ei +E; =0. (8-2) 


‘Combination of these two equations yields 
(/forr.R9i,@R) aR’) = -£/(R). (83) 
' 
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This is the fundamental equation that we wish to solve with the 
method of moments. 

First, the Green’s function must be derived. The Green’s func. 
tion derived here is for scattering in free space. The waveguide iris 
discussed in Chapter 9 offers another example of scattering, with 
the Green’s function changed by the waveguide boundary condi- 
tions. The scattered electric field can be expressed in terms of the 
magnetic vector potential as 


E‘= [V(v - A) + B74], (8-4) 


JOM Eo 


where in the usual notation 6? = ye. In turn, the magnetic 
vector potential is given by 


}D 
A(R) = [f= Hol, sa em @R, (85) 


where the distance D is 


D=|R-R|=V(x-xP + (y-yP+(2-zF (66) 
in rectangular coordinates. Combination of (8-4) and (8-5) yields 


e IBD —ipD 


-- aamal| [r- LE) >| + BUR) 


aR. 
(8-7) 


The various forms of this basic equation can be solved in the usual 
way with the moment method [1]. In contrast to most of the 
relatively simple examples we have considered in earlier chapters, 
however, the integrals that are be evaluated to define the matrix 
coefficients must be evaluated numerically. The result is that for 
scattering and other difficult problems, the majority of the com- 
puter run time is devoted to “filling the matrix,” that is, calculation 
of the matrix terms. 
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8.2. SCATTERING OF A PLANE WAVE OFF A THIN CYLINDER 


For a specific example, consider the uniform plane wave traveling 
jn the x-direction: 

Ei(x) = £,,e7*. (8-8) 
‘This wave encounters a thin cylinder (wire) of radius a and length 
L, whose axis lies on the z-axis. If the wire is sufficiently thin, then 
approximately on the surface of the cylinder, 
EX(x) = E, (8-9) 
For this situation, the current flows only in the z-direction and is 
uniform around the cylinder, and hence has the same effect as a 
current concentrated on the z-axis. Thus the z-component of the 
vector potential can be written in terms of the total current as 


Ho 


Afr, j-7f 5 m2!) de’, (8-10) 
where in the cylindrical coordinates used in this equation 

D=yr?+(z-2'). (8-11) 
On the surface of the cylinder 

Ho Aid 
AJa, Ded ie pie") ae", (8-12) 

where now 

D=yat+(z-2'). (8-13) 


Because D depends only upon z — z’, it is convenient to write 
e7iBD 


D 


K(z-2')= (8-14) 
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and (8-12) can be written as where now the Green’s function is 


L & 2 ’ 
A a,z)= mand K(z —2')1(z') dz’. (8-15) G(z,2')= (5 +e) [K 2) + K(z+z’)]. (8-22) 
Thus, on the surface of the cylinder, (8-2) and (8-4) reduce to With this Green’s function, the integral equation (8-17) can be 
solved with the moment method. The integration required to obtain 
the matrix coefficients must be performed numerically. Useful re- 
sults can be obtained using pulse functions for the basis functions, 
but as the number N of basis functions is increased, these integrals 
do not converge well and the behavior of the solution as a function 
of N is not too smooth. Basis functions that correspond to a 
smoother approximation for the current produce better results. 
Another approach can be taken to improve the convergence of 
the solution process. Equation (8-16) is a second-order differential 
equation, and thus the integral can be written as 


F} # 2 eb , , , 
—j4mwe Ey = (5 +6 | I K(z-z')K(z') dz’. (8-16) 


Several approaches have been used with this equation. The direct 
one is to apply the operator to the integrand and get 


Sit metp ee [o, 2')I(z') dz', (8-17) 


where the Green’s function is fe 1 
4 f G(z, 2’)I(z') dz’ = PAT eck m + Bsin( Bz) + Ccos(Bz), 
a 0 

G(z,2') = (= eS a -z'), (8-18) (8-23) 
which is one form of Hallén’s equation. This integral equation can 
be solved with better convergence properties than Pocklington’s 
equation. The coefficient B is seen to be zero because of the 


which is one form of Pocklington’s equation. Performing the dif- 
ferentiation leads to the explicit result 


2iB 2+a'B? 3jBa? 3a? symmetry, and C must be determined when the current is calcu- 
G(z, 2’) =e?) —> + SS ar > ers (OD lated. If N —1 basis functions are used and N conditions are 
D D D D satisfied, then solution of the system of N equations yields the 


N-—1 parameters that define the approximation to the current. 
One can make use of the fact that the current is zero at the end of 
the wire. 

Butler and Wilton compare in [2] various numerical techniques 
that have been used with this scattering problem. 


The current can be seen to be symmetrical around the midpoint 
of the wire, and thus if the origin is redefined so that the equation 
(8-17) becomes 


=j4rweyE,, = pa 7G Z)I(2’) az', (8-20) 


8.3. SCATTERING OF A PLANE WAVE OFF A THICK CYLINDER 
the fact that J is an even function results in 
Tf the cylinder is thick, which means that the radius a is not small 
compared to the wavelength of the incident wave, several of the 


L/2 
Aw we En = [© G(z, 2')I(2') a", (8-21) red te : : n al of | 
0 approximations made in the preceding section are not valid. First, 
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(8-9) is not a good approximation. If we write the incident field in 
cylindrical coordinates (r, @, z), then 


x=rcos@ 
y=rsind, 


(8-24) 


and the incident field can be written from (8-8) as 


El(r, 8) = Eger? (8:25) 


Second, the current depends upon the angle # around the axis and 
(8-10) becomes 


aby 


Ar, $,z) = = 


[i [remy (a ¢ z') dd! dz' (8-26) 
0-0 pee y 


where 


D = yr? — 2ar cos(s — 6!) +a? + (2-2). (8-27) 


The electric field corresponding to this vector potential is, as 
before, 


1 a 
2 


Joe ty a +B) Ar dy2) (8-28) . 


Ex(r,b, 2) = aa 


Combination of these equations and setting r = a gives the desired 
integral equation, 


a 


—E,,e7 184005 4 


ia 


1 
X HIG. 2!) db dz’, 


J4mwey 


(8-29) 


J 
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which corresponds to (8-19) and where now 


(8-30) 


In addition, the current on the ends of the cylinder, which was 
previously neglected, now must be included. 


PROBLEMS 


1. Use integration by parts on equation (8-17) to derive another 
form of Pocklington’s equation 


jae Ee = [ke —2')[1"(2') + B(z')] de’ 
4 AK(z) + BK(z -L) 


where A and B are constants determined by the current J. What 
are A and B? 


2. Using the fact that B = 27/A and the change of variables 
a@ = Bz to show that the nature of the solution to Pocklington’s 
equation depends upon only the two parameters a/A and L/A. 


3. Derive the integral equation corresponding to (8-20) if the inci- 
dent wave is not normal to the wire. 


COMPUTER PROJECT 8-1 


Write a computer program to solve with the moment method the 
integral equation described by equations (8-21) and (8-22). Use 
Pulse functions for basis functions and point match at the center of 
the pulses. Run the program for the following parameters: 
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and plot the solution for the current at the center of the wire as 
function of 1/N, where N is the number of basis functions. Do i 
solutions converge well as N is increased? 7 


REFERENCES 


1. R. C, Hansen (ed.), Moment Methods in Antennas i 
Artech House, Norwood, Mass., 1990. ii 
2, Chalmers M. Butler and D. R. Wilton, “Analysi i 
n ; |D. R. a ysis of various numeri 
techniques applied to thin-wire scatterers,” JEEE Transactions one 
tennas and Propagation, AP-23, 1975. 
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9.1. BASIC FOURIER SERIES RELATIONS 


In practice the most difficult task in the application of the method 

moments usually is determination of the Green’s function. Only 
in a few situations, such as the parallel-plate capacitor in un- 
bounded space, can analytic expressions for the Green's function be 
derived. In principle, as shown in Chapter 7, determination of the 
Green’s function is straightforward. If the function to be found is 
‘the solution of an operator equation, as in 


[f(R)] = 8(R), (9-1) 
with the boundary condition 
f(R) = 0 on the boundary, (9-2) 
the Green's function G(R, R’) is the solution of 


with the same boundary condition. 


Z[G(R,R’)] = 3(R - RF’) 
(9-3) 


Sometimes, the defining equation (9-3) can be solved directly in 
closed analytic form. For many situations, however, solutions for 
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the Green's function as an infinite series or as an integral transform 
are required. Waveguide and shielded microstrip examples are 
good illustrations of such a situation and will be considered later in 
the chapter. 


9.2, EXAMPLE INVOLVING LAPLACE’S EQUATION 


First, consider a simple example, with a variable E(x, y) that 
satisfies Laplace’s equation inside the rectangular region shown in 
Fig. 9-1. The function is zero on three sides and equal to f(x) on 
the upper side. If for each value of y the function is expanded in a 
Fourier series, 


N 
E(x,y)= 0 A,y)sin =. (9-4) 


n=l 


If each term in this series satisfies Laplace’s equation, 


(9-5) 


FIGURE 9-1. A two-dimensional potential region. 
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To satisfy the boundary conditions, the appropriate solution is 


‘ nwy 
Ay) = @, sih—, (9-6) 
and the series is 
us nvy | nwx 
E(x,y)= Y a, sinh — sin——. (9-7) 


n=l 


1 series satisfies the zero boundary condition on three of the 
les. To satisfy the boundary condition on the upper side, the 
oefficients a, are determined in the usual Fourier series manner 


ge 2 .@ nix’ ") ax’ 0.8 
a, sin oo ay sin f(x (9-8) 
the series is 
N sinh(nwy/a) . nx 2 .a. nix’ 
y= eH ") dx’. (9-9 
ey) » sinh(nzb/a) a ais a F) 9) 
expression can be written as 
a 
B(x. y) = [ G(xy,x')fl2') ae’, (9-10) 
the Green’s function is 
2 % sinh x j 
Baancimee Ee sinh(nay/a) | 7. ant | (6-11) 


Ce sinh(nzb/a) a 


the boundary conditions are nonzero on the other three sides, 
three other similar series must be added to the solution. 
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9.3. INDUCTIVE IRIS IN A WAVEGUIDE 


An example of a scattering problem that is also a good illustration 
of the use of Fourier series techniques is the inductive iris in a 
waveguide. With the usual notation for a rectangular waveguide, 
the TE, mode has an E-field defined by ; 


nwx 


E,,(%, 2) = sin I ee, (9-12) 
where the propagation constant y,, satisfies 
4 aE P 2 2 2, 
% = (= =e5 where k* = w*e. (9-13) 
Usual waveguide practice results in 
T 2a 
ee (9-14) 


and thus y, is imaginary and y, is real for all values of n greater 
than 1, Now assume that a dominant mode wave defined by 


o WX 
E\(x,2) = sine" (9-15) 


is traveling in the positive z-direction and at z = 0 there is a thin 
diaphragm (iris) with the cross section shown in Fig. 9-2. If this iris 
is a conductor, current will flow in the y-direction and create a 


0 c a x 


FIGURE 9-2. Inductive iris in a rectangular waveguide. 
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ttered wave. If the iris is a perfect conductor, the total E-field 
the iris will be zero, To determine the scattered field, we need to 
ind the current on the iris. First we need an equation that relates 
1e scattered field to the current. Although the incident E-field is a 
nusoidal function of x, the resulting current is not, The current 
density can be expanded in a Fourier series of the form 


Rel Sea, (9-16) 
n=) @ 
where the coefficients c,, are given by 
2 ,@ nx’ 
~ i sin——J,(x') de’. (9-17) 
‘or z slightly less than 0, the magnetic field is given by 
Hi(x,0—) = —,(x). (9-18) 
ause for each mode the E and H fields are related by 
1, (9-19) 


‘n 


hd al nx 
E;(x,0) = (- sin-—— 9-20 
vl > 2, a ( ) 
for z = 0, and for z <0 
Be 2)= 2b (9-21) 


n=l 


144 
Substitution of (9-17) and rearrangement yields 


© jon | nwx 


EX(x,z)= — sin ene aig ef x’) dx’. 
i adil A Dr 2’) 
In order that the total field on the iris be zero, 
2 _ 7x 
E}(x,0) = ~sin=—. 


which leads to the integral equation 
© 7x 
i G(x, x')J,(x") de! = sin—- for0 <x <c, 
0 


where the Green’s function is 


, 


= Jom  nwx . nx’ 
Cz2,2)= ¥ sin—sin 
n= 4¥n a 


The moment-method solution based on the approximation 
N 

A(x) = L aihj(x) 
jet 


and the weighting functions w,(x) leads to 


or aeane ¥ ies ie 3 _ 7x 
Kal, [miele x W(x) ae de = f w(x)sin—— 
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dx, (9-27) 
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For the Galerkin method with pulse functions, use N subintervals 
of length 


e 


9- 
(9-22) ee 


h= (9-29) 


and the first integral is 


nTxX 


(9-23) [wdesin a =f" ae (9-30) 


which easily can be evaluated as 


a [ox = »| = cos“ 


(9-31) 


(9-24) ia w, (x)sin dx = 


This integral depends on n and j and the parameter h/a, which 
can be expressed as 


(9-25) 
(9-32) 


9.26 The other two integrals are of the same form and thus it is 
(9-28) convenient to write (9-28) as 
Nv = jw 
> ay ys 


j=l n= Yn 


= ( \'F(n, iF(n,j) = <r, i). (9-33) 


. nzh nzh | 
which can be written F(n,i) = cos “= - 0| - cos ). (9-34) 
a “ jap n7Xx <. nwx’ ined 5 . 5 
x a0 = Uf w,(x)sin = as| Lf sin— W(x’) de Tt is convenient to rearrange (9-33) as 
c wx N wma = Gis 
= [ w)x)sin— de (9-28) bale (1,i), (9-33) 
0 a 
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and then by introducing the normalized set of parameters which is 
_ ema N jou ¢, 7x 
X,= a (9-34) p=- care, a A) dx. (9-38) 
the equation becomes In terms of our previous notation, 
z Sf z ExF( i) (9-39) 
X, DL —F(n, i)F(n, j) = F(, i). - Oe ee 
i=l Py nay, (FC J) = F144) (9:33) ay) j=1 


In circuit terms, the effect of the iris on the dominant mode can be 
represented as a shunt impedance across a transmission line. The 
normalized shunt impedance can be shown to be 


The coefficient ay,, is found from (9-13) as 


ay, = V (nr) — (ka). (9-34) 


ZL l+p 

i ———, (9-40) 
Because for large n this coefficient is proportional to only the first Zo 28 
power of n, the series for the Green’s function can be seen from 
(9-25) not to converge well. The two integrals, however, each The right side of this equation can be shown to be an imaginary 
introduce a factor of 1/n, and therefore the terms in the series in number, leading to the iris being referred to as an inductive iris. 
(9-33) go to zero as 1/n* and therefore the series converges well. 

The dominant term in the scattered (reflected) wave is the n = 1 

term in (9-22) and thus the reflection coefficient at z = 0 is 9.4. BASIC FOURIER INTEGRAL RELATIONS 


jou px In Section 9.3 we demonstrated the usefulness of Fourier series 

p=-— f sin—J,(x) dx. (9-35) representations for variables defined on finite domains, as in the 

ayy “o b case of waveguides and shielded microstrip circuits. For variables 

¢ defined on infinite domains, the series are not useful and we must 

Using the basis function expansion we used, turn to the Fourier transform to use similar methods. The Fourier 
transform F(q) of a function F(x) is defined by the integral 


N 
Fu = 1-36 o + 
y(e) LaF). (33g) F(a) = [eM F(x) de. (9-41) 
the reflection coefficient can be expressed as The function F(x) can be recovered from its transform F(a) by the 
integral 
jon -¢. wx NX 
= —— | sin— 9-37 1 ,+= bs 
p rae Rs LeasPi(x) de, (9-37) F(x) = xf eveF (a) da. (9-42) 
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These relations are analogous to the Fourier series. The function Taking the inverse transform gives 
F(x) now is defined over an infinite domain rather than the finite +s 
domain in the case of the series. The function F of the continuous E(x,y) = f G(x, y, x’) f(x") dx’, (9-46) 


variable @ replaces the discrete set of Fourier coefficients in the 


Fourier series case. . - 
where the Green’s function is 


ith _,) Sinh ay ay 
9.5. FOURIER INTEGRAL SOLUTION G(x, yx!) = 5 ale Smbop**™ = (9-47) 
Consider now the situation where the side walls in the structure The integral in this expression probably must be evaluated numeri- 
shown in Fig. 9-1 are removed to leave the region shown in Fig. 9-3, cally. 
Suppose that the function E(x, y) is zero for y = 0 and is equal to 
f(x) for y = b. The transform of both sides of Laplace’s equation is 
PROBLEM 
ah ss 
[-« + }é(a y)=0, (9-43) 1. For the inductive iris, derive expressions for A,, and B, for the 


basis and weighting functions 


which has solutions that can be expressed as linear combinations of _ iwx 
hyperbolic sines and cosines. Because of the first boundary condi- U(x) = ee 
tion, the sinh function will be used and thus ar 
4 w,(x) = sin—. 
E(a, y) = ¢(a)sinh ay. (9-44) ts 


Use of the second boundary condition to evaluate c(@) leads to 


‘ COMPUTER PROJECT 9-1 
sinh ay 


E(a,y) = F(a). (9-45) 


sinh ab 


Write a program to implement the Galerkin with pulse functions 
Solution of the inductive iris problem. Use the value for ka of 4.5, 
and solve for various values of c/a and graph p as a function of 


c/a from () to 1. 
26, te _ REFERENCES 
2 
an?” ay? { 1. D. Mirshekar-Syahkal, Spectral Domain Method for Microwave Inte 
— grated Circuits, John Wiley & Sons, New York, 1990. 
O o=0 


: 2. Craig Scott, The Spectral Domain Method in Electromagnetics, Artech 
FIGURE 9-3. A two-dimensional potential region without sides. House, Norwood, Mass., 1989. 


CHAPTER TEN 
mee 


Spectral Analysis of Microstrip 
Transmission Lines 


10.1. MICROSTRIP TRANSMISSION LINE 


o versions of microstrip transmission line will be analyzed here, 
lded and unshielded, shown in Figs. 10-1 and 10-2, respectively. 
scause of the two valucs of permittivity, microstrip transmission 
es cannot support a TEM mode of transmission. At sufficiently 
frequencies, however, the dominant mode is approximately 
, and quasistatic or quasi-TEM calculations based on the 
ady-state values of inductance L and capacitance C per unit 
length give useful results. Most analyses do not calculate L directly. 
nstead, the capacitance C, that results from assuming that €, = €) 


e true value of €,. These two values of capacitance are then used 
determine the transmission line parameters such as phase veloc- 
ity. As the frequency is increased, the quasi-TEM results become 
less accurate, and eventually the “full-wave” solution that results 
using Maxwell's equations must be determined. 


10.2, QUASI-TEM SPECTRAL ANALYSIS OF SHIELDED 
ICROSTRIP TRANSMISSION LINE 


Consider now the quasi-TEM solution for the shielded microstrip 
transmission line. The capacitance per unit length is to be calcu- 
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and 


#(0.5a, y) = 0. (10-3) 


_ The potential in the lower region (y < h), which will be denoted by 
@,, also satisfies the boundary condition 


,(x,0) = 0. (10-4) 


The separation of variables approach leads to the solution 


n nwTxX 
(x,y) = Lie, sinh a cos =, (10-5) 


za 


~ 


_ where the prime indicates summation over odd positive values of n 
(1, 3,5,...). The potential in the upper region (h < y), which will 
_ be denoted by #p, also satisfies the boundary condition 


FIGURE 10-1. Shielded microstrip transmission line. 


lated. In regions where e is constant, the potential satisfies Laplace’s 
equation 


{ (x, b) = 0, (10-6) 


and the separation of variables approach leads to the solution 


rb Pb 0 (10-1) ee ae 
Pee re ag bs 
ax? * ay? (x,y) = D'd, sinh "cos “= (10-7) 


The charge density, which is at the interface, can be expanded in a 


The potential clearly is an even function of x and on the side walls 
Fourier series 


satisfies the boundary conditions 


NTX 


$(—0.5a,y) =0 (10-2) p(x) = Lf, cos ——, (10-8) 


where the coefficients can be written as 
f 2 | om niax a (10-9) 
=— x)cos —— dx, a 
ay ee a a aa 


The coefficients ¢, and d,, will be determined by satisfying term by 
term the conditions at the interface between the regions (y = h), 
The continuity of the potential 


$o(x-h) = b\(x,h) (10-10) 


4 


FIGURE 10-2. Unshielded microstrip transmission fine. 
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leads to 


nth nw(h —b 
e, sinh ~~ =d, ith 


and the interface condition 
a a 
£055 Po - S151 = Py 
at y =h leads to 


nt nw nw nr 
ed cosh = th —b)+ en cosh =k Fhe 


Combination of these equations leads to 


+050 nax' 


nax 
=Y' — ——pi(#') a's 
H(x,h) 8, 605" f cos ——p,(*") ¢ 


—0Sa 
where 


2/nwr 
—e, coth(na/a)(h — b) + €, coth(nm/a)h” 


8, = 
which can be written as 


(ash) = [7 G(x, xox’) de 


—0. 


where the Green’s function is 


nwx 


G(x, x’) = Dg, cos = 


If the upper shield is removed, or alternatively, taking b = ~. the 
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term g, becomes 
— — 10-18 
(10-11) Bn e+e, coth(na/ayh pacts) 
If the potential on the strip is taken to be unity, (10-16) gives the 
integral equation that determines the charge density as 
(10-12) +0.5a 
J G(x, x')p,(x') de’ = 1. (10-19) 
~OSa 
10-13 10.3. QUASI-TEM SPECTRAL ANALYSIS OF UNSHIELDED 
(10-13) MICROSTRIP TRANSMISSION LINE 
Tn terms of the Fourier transform 
- +2 
(10-14) Hay) =f eM*b(x, y) de, (10-20) 
Laplace’s equation becomes 
a = 
(10-15) | 57 «}é(ay) =0. (10-21) 
dy? 
The solutions to this equation can be expressed in terms of expo- 
nential functions or in terms of hyperbolic functions. To satisfy the 
(10-16) boundary condition that the potential is zero at y = 0, the solution 
for y <h is 
b\(a, y) =c(@) sinh ay, (10-22) 
(10-17) and for h < y the solution is 
do(a, y) =d(a)e-*”. (10-23) 
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Continuity of the potential at y = h gives 


d(a)e"s* 
ath) = sinhah * (10-24) 
and the interface condition 
a a 
#05, Po at “15%! = Py (10-25) 
at y =/ leads to 
—e,d(a)a@ + €,c(@)a cosh ah = p,. (10-26) 
Combination of these equations leads to 
+2 
(xh) = f G(x, x’)p,(x’) de’, (10-27) 
where 
1 ,+# 1 
"= a2) da. (10-28 
Slax) ro a(ey + €, coth ah) em ) 
Use of the symmetry with respect to @ gives 
G(x,2") = + [cos hie ——>——_ ja. Wael 
, ay a(e) + €, coth ah) 


10.4. COMPARISON OF FOURIER SERIES AND FOURIER 
INTEGRAL SOLUTIONS 


Physically, one would expect that if the shield over the shielded 
microstrip were far from the strip that the solution would be 
approximately the same as the solution of the unshielded mi- 
crostrip. Mathematically, this means that the Fourier series and the 


COMPARISON OF FOURIER SOLUTIONS: 157 


Fourier integral solutions would be related. One way to see this is 
to examine the Green’s function for the shielded microstrip. 


G(x, x') = = [cos a(x —x') de. (10-30) 


a(e, + €, coth ah) 


The integral can be evaluated numerically by dividing the region 
from 0 to = into subintervals of length 26 and using the midpoint 
integration formula. The points of the intervals are nd, where n is 
an odd positive integer. This evaluation of the integral results in 


1 1 
G(x,x')=—)’ 5(x — x’) ———_——_—____26. 
(eee) 7 Psemnile *) ae, + €,coth n 5h) ie 
(10-31) 
Comparison with the Fourier series solution 
nix nix 
G(x, x') = D'g,, cos ag (10-32) 
where 
2/nizr 
= —— (10-33) 


€) + €, coth(naz/a)h * 


shows that the solutions almost agree if the subinterval in the 
numerical evaluation of the integral satisfies 


5- (10-34) 


a|a 


If the width a is large, the interval 28 is small enough that the 
Numerical evaluation of the integral should be accurate. One dif- 
ference remains because the trigonometric term in the integral is 
cos a(x — x'), whereas the corresponding term is cos ax cos ax’. 
Clearly, in the limit the term involving sin ax sin ax’ must sum to 
zero. 
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PROBLEM 


1. Derive the result in Section 10.2 for no top cover (b = =), 
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CHAPTER ELEVEN 
FS 


Spectral Analysis 
of Microstrip Circuits 


11.1. MICROSTRIP CIRCUITS 


In this chapter we consider the spectral analysis of microstrip 
circuits. The basic structure is similar to that of the microstrip 
transmission line considered in Chapter 6. The substrate sits on a 
ground plane, but the metallization on the top of the substrate no 
longer is restricted to a constant-width line but can have an arbi- 
trary shape in the x —z plane. The methods used are similar to 
those of Chapter 10, with the difference that the lateral variation is 
in two dimensions rather than one. As in Chapter 10, we consider 
both shielded and unshielded microstrip. We first perform qua- 
sistatic analysis for both types of circuit and then do full-wave 
analysis for the unshielded microstrip. The full-wave analysis of the 
shielded microstrip circuit is similar except that Fourier series 
rather than Fourier integrals are used. The quasi-static analysis 
calculates the capacitance with respect to the ground plane of an 
arbitrarily shaped metal figure on the substrate, whereas the full- 
wave analysis determines the relation between current in the mi- 
crostrip and the electric field across the substrate. 

Problems and computer projects are not included at the end of 
this chapter, because writing a computer program to implement 
these methods is a major effort and the results are too close to the 
level of current research in this field. 
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11,2. QUASI-TEM SPECTRAL ANALYSIS OF SHIELDED 
MICROSTRIP CIRCUITS 


Consider now the quasi-TEM solution for the shielded microstrip 
circuit, which is a three-dimensional problem. In regions where e is 
constant, the potential satisfies Laplace’s equation 


76 Fh Hh 
ate * ge (11-1) 


The separation of variables approach leads to the solution in the 
lower region of 


nwz 


nTx 
(x, ¥.2) = Yen, Sinh(u,,,¥)sin aa sin ee (11-2) 
mn 
where 
mm \? nt \? 
nn = (=) +(=) > (11-3) 
a c 


The potential in the upper region (A < y), which will be denoted by 
by, can be written 


naTXx ATrzZ 
bo(X. 9.2) = Ld, sinh(u,,,(¥ — b))sin— — sin =e (11-4) 


mn 


The charge density, which is at the interface, can be expanded in a 
Fourier series 


mrx nwz 
p(x,z) = L finn sin sin) (11-5) 


mn a 


where the coefficients can be written as 


4 pea max nz 
=— i ‘in — , 11-6 
‘a aed, fea 2)sin 7 i dxdz. (11-6) 
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The coefficients c,,,, and d,,,, Will be determined by satisfying term 
by term the conditions at the interface between the regions (y = h). 
The continuity of the potential 


by(x.h, 2) = (x, A, z) (11-7) 
leads to 
Cm Sinh(u,,,, ¥) = d, Sinh(u,,,(¥ — b)), (11-8) 
and the interface condition 


a a 
£05, % + e159 =p, (11-9) 


leads to 


= € 4 nnn COSD teyyy( ft — B) + € 1Cyl mn COSH Mina = finns 
(11-10) 


Combination of these equations leads to 
c pit 
(xh, 2) = ff G(x, 2,2", 2')p(2", 2") de’ de’, (11-11) 
0-0 


where the Green’s function is 


ra Sak yr _ max  nwz  onwx’ — nz’ 
.2/y= — sin —— sin —— 
(er an | Pe Sen sim z sin z sin — 3 
(11-12) 
and where 
T/tyn)(4/ac 
es (1/tinn)(4/a0) tie 


—e, coth u,,,(h — b) + €, cothu,,,A° 


The solution using the moment method then proceeds in the usual 
manner. 
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11.3. QUASI-TEM SPECTRAL ANALYSIS OF UNSHIELDED 
MICROSTRIP CIRCUITS 


Now consider the unshielded microstrip circuit. As before, at low 
frequencies the electric potential is governed by Laplace’s equation, 


rh Fb #d 


aa? * ay * age ite 


In terms of the Fourier transform of ¢, defined by 
H(a,y.B) = [fee h(x, y,z) dedz, (11-15) 


where unless otherwise indicated limits of integration are —* and 
oc, Laplace’s equation becomes 


# a 
[-«+ ae -#°}6 ~0, (11-16) 
which is 
# 
gr-# bo, (11-17) 
where 


u=yar+ B. (11-18) 
In the dielectric material, where y < h, the appropriate solution is 
$,(a, vy, B) = ea, B)sinh uy (11-19) 

and above the interface 


do(a.y,B) = d(a, Be. (11-20) 
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Because the potential is continuous across the interface at y=h 
and therefore its transform also is 


c(a@, B)sinh uh = d(a, B)e~"". (11-21) 


Because the interface condition on the normal D leads to 
i) Cl) 
£09, Fo + e15y% =f (11-22) 


the transform satisfies 


B x @ = 
~ egy bo + 615541 = fy (11-23) 


Substitution of the two expressions for & gives 


egude~“* + €,uc cosh uh = p,. (11-24) 


Combination of these equations gives for the transform potential at 
the interface 


= 1 
aH) a othe wet (11-25) 


The potential @ can then be found as the inverse transform 


1 


hep heya aay eB Se dp. 


(11-26) 


F(z, hi, 2) = (ys logit 


Substitution of 


Aya. B) = ffe-*e*p,( x’, 2") dx" de" (11-27) 
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and interchange of the order of integration yields 
Px, h, z) = ikaes 2,4", z')p,(x", 2) de’ dz’, (11-28) 


where the Green’s function is 


1)? sic teeecagl 
G(x, z,x',2') = (=) [peme Ye IB2-2') 


Tee coth any W- (11-29) 


With the variable changes 


@=ucosw 
B=usiny 
x —x'=rcosé 
z-2’=rsin@ (11-30) 


the Green’s function becomes 


1 


2 0 2. 
G(x, 2, 2,2’) = (=<) f i, ee 
0-0 


See . (11-31 
* Te, + €, coth 7 aia ( ) 


Because the Bessel function of zero order can be represented by 
the integral 


1 pan. 
= — f“eireosegg, 11-32) 
42) = z= fe ( 
the integral with respect to 6 is 


[Cem do = 2aJy(ru), (11-33) 
0 
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and the expression for the Green’s function becomes 


1 p= Jo(ru) 
ey) el €) + €,coth uh ae (a4) 


Because the variable r is given by 


r=V(x-x'f t(z-2'V, (11-35) 


the Green’s function depends only on the distance between the 
primed and unprimed points. The integral equation can be solved 
with the moment method in the usual way, but numerical integra- 
tion clearly will be required to determine the Green's function. A 
useful approach is to evaluate and store values for this function at a 
discrete set of values of A and when the function is needed for 
other values to interpolate between the discrete values. 


11.4. FULL-WAVE SOLUTION OF UNSHIELDED 
MICROWAVE CIRCUITS 


We shall derive the integral equation by first using the fact that the 
electric field satisfies the wave equation. In a region where x and € 
are constant, the E-field satisfies 


5 +k)E=0, (11-36) 


where as usual 
k? = wpe. (11-37) 


The two-dimensional Fourier transform of E, defined by 


E(a, y, B) = f fete E(x, y, 2) dedz, (11-38) 
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satisfies 


td 
[-« tee ~ p+) =0, (11-39) 


This equation can be written as 


aro |e =o, (11-40) 


where 


u= ya" +p? ~k?. (11-41) 


This is a second-order ordinary differential equation with respect to 
y and has solutions that can be expressed either in terms of 
exponential functions or in terms of hyperbolic sine and cosine 
functions. In the lower region where y </, the x and z compo- 
nents are equal to zero on the ground plane and thus are propor- 
tional to the hyperbolic sine function. To simplify later calculations, 
these components can be expressed as 


= ___ sinh(w,y) 
E,(a, y,B) = Seah (ash) (11-42) 
and 
E,(a, y,B) = ce, (11-43) 
1 


where the ¢ coefficients are functions of a and B. Because the 
divergence of the electric field is zero, the y component can be 
evaluated as 


Jac, + jBc, cosh(u,y) 


————, (11-44) 
uy sinh(u,h) 


E,(a, yB) = 
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In the region above the interface (y > h), the transforms of the x 
and z components of the electric field can be written 


Eq = cet (11-45) 


and 
E.=c,e""!™, (11-46) 


where we have used the continuity of the tangential electric field 
across the interface. Again, the fact that the divergence of the 
electric field is zero can be used to show that 


2 jac, + jBc, 
>= 


uy) 4 
ie € (11-47) 


From the Maxwell curl-E equation, the transforms of the tangential 
magnetic field components can be shown to satisfy 


‘ Oo e 
jou, = ~ > E, + iBE, (11-48) 
and 
5 4 ao x 
jopH, = —jaE, + yet (11-49) 


The components below the interface can then be evaluated as 


aBe, + (uz — B*)c, cosh(u,y) 
i: uy sinh(u,h) 


jopH,, = (11-50) 


and 


out, = WE Te }ee = @Be, cosh(uy) 
eS uy sinh(u,h) ” 


(11-51) 


and the corresponding components for the magnetic field above the 


168 SPECTRAL ANALYSIS OF MICROSTRIP CIRCUITS 


interface are 
i —aBe,(uji — B)c 
jopAyy = abe. (ub PN ys (11-52) 
ug 
and 
2 uz — a*)c, — aBe, 
jopH. = WHE Vee BPE aya (11-53) 


uo 
At the interface, as y approaches A, these components simplify to 


ses) 2_ g2 
jap j(avh,B) = SS ES tii (11-54) 
1 


(uz —a*)c, — aB 


jopH.(a,h,B) = = coth(u,h), (11-55) 


uy 
= 2_ g2 
jopl,o(a,h, B) - Ses, (11-56) 
0 
and 
=: _ 
pti sy= ae (11:57) 


ug 


The tangential H-field is discontinuous at the y =f interface in 
accordance with the equations 


(11-58) 
and 
-Hyyt Hy =F, (11-59) 
which lead to 
(uz — a?)e, — ae. A (uj —@7)e, — aBe, slit, 
Uo by 
=—jopd, (11-60) 
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and 
—afe, + (u2 — B*)c. tg (u? — a)c, — ape, weit 
uy uy 
= -jopJ.. (11-61) 


These are two linear algebraic equations in the two unknown 
parameters c, and c., and making use of the relations 


¢, = E,(a,h, B) (11-62) 


c, = E,(a,h, B), (11-63) 
the results are 
—jwE,(a,h, B) 


u, + 4, coth uh 
€o(U, + Ue, coth ujh) (uy + u, coth u,h) 


- Me 
+ @ J, ————_—_ 
Ug + u, coth uh 


= (a7, + aBJ.) 
(11-64) 


and a similar equation for E_. With the notation 


é u, + Ug coth uh 11:65 
© g(t, + Upe, coth uy) (uy + w, coth ujh) ( ) 


o m 
i 11-66 

"ty + u, cothu,h ( ) 
the equations can be written 


—jwE,(a,h,B) = (aS, + oBJ.)G(a, B) + w°S,G,,(a, B) 
(11-67) 
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and 
~jwE(a,h,B) = (Bs, + B*E)G(a,B) + w°S.G,(a,B). 
(11-68) 
In taking the inverse Fourier transform, the factors a and B 


introduce derivatives, which can be handled in several ways. One 
way is associate the derivatives with the entire term, which leads to 


a 
—jwE,(x,h,z) = — pe [J Galr eC, 2") a dz’ 
fg 
~ ax dz 


of wo [G(r J," 2’) (11-69) 


[ Gr) _(x, 27) de! de! 


and a similar equation for E.. A better way numerically probably is 
to associate at least one of the derivatives with the currents J, and 
J., When the moment method is used, this requires the use of 
continuous basis functions. If both derivatives are associated with 
the currents, the basis functions must have continuous derivatives. 
With one derivative associated with the currents, 


a ad, 
-jwE,(x,h,2z) = ~ A SfEMz 2') dx’ dz! 
a as, 
— ay [JG G A" 2) de! ae! 
+0? [[G,,(r)J,(2", 2") dx’ do’, (11-70) 


This equation and the corresponding equation for E. can then be 
written in vector form as 


—jwE(x,h,z) = Vf [GV + J(x', 2) dx’ dz’ 


+ 0? f [Gy(r)I(x’, 2") dx’ dz’. (11-71) 
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These equations frequently are written in a different form. The 
form selected here is consistent with 


= —Vo —jwA, (11-72) 
which in the x-direction is 
ewe, 
PA aS —jwA,. (11-73) 


The Fourier transforms of these components are 
E, = -jad -jwA,, (11-74) 


where the transforms of the potentials are 


6 =G.p, (11-75) 
and 
A=6G,J, (11-76) 
which leads to 
E, = ~jaG.p, — joG,J, (11-77) 
because 
: a, ad, 
i ala ae ae (11-78) 


which in terms of the transforms is 


~ieS, — ips. 


jo 


p= (11-79) 
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Substitution of this expression gives 


= 6 os, a aps, 


E,= GF — jd. (11-80) 
Comparison with 
Bes. § al, + aps. 
sles) = — 
u, + uy coth uh 
x 
€y(u, + Uge, coth uyh)(uy + u, coth uh) 
2 B 
is 11-81 
hil uy + u, coth uh ( ) 
shows that 
P 1 uy + ug coth uh 
@. =. eee (1189) 
. jw e (u, + ue, coth u,h)(uy + u, coth wih) 
and 
3 Me 
8 Set 11-83 
Gm Uy + uy coth uh ( ) 


CHAPTER TWELVE 
De 


Mode Matching 


12.1. APPROXIMATING A FUNCTION BY A SERIES 


A fundamental problem in matching is approximation of a function 
that is approximated by a finite series, which in one dimension can 
be written 


N 
g(x) = Yo a)B(x). (12-1) 
j=0 


Frequently, the finite series is obtained by truncation of an infinite 
Series that should equal the function. The problem is to select the 
coefficients g, so as make the approximation in some sense a good 
approximation. We have already encountered this problem as part 
of the moment method, but it is encountered in a wider context. A 
fairly general approach is to select a set of N +1 weighting 
functions w,(x), multiply both sides, and perform an integration 
over the domain of x, resulting in 


N 
Joule(x) de = ¥ 8, fw,(x)B,(x) de. (12-2) 
j=0 


As discussed in Chapter 7, this method is frequently referred to as 
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the method of weighted residuals, since if the residual is defined as 
N 
R(x) = (x) — Li g:B(x), (12:3) 
i=0 
the weighted residual is equal to zero if 
N 
J mo] (2) =x 28,0] dx = 0, (12-4) 
j=0 


which is equivalent to the equation derived earlier. A specific set of 
weighting functions sometimes used is a set of delta functions 


w(x) = (x —x,), (12-5) 


which leads to the point-matching result 


N 
a(x;)= Y 3/B(x,). (12-6) 


j-0 


The Galerkin selection of weighting function was another approach 
discussed when the moment method was covered. A similar ap- 
proach is the least square error (LSE), where the coefficients are so 
selected as to minimize the integrated squared error 


2 
ISE = je = E 5,84) dx. (12-7) 
j=0 


The minimum value of this expression is defined by the set of 
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equations 
a N id 
wa Ss - Lenco] dx = 0, (12-8) 
which leads to 
N 
[Blx)a(x) de = D 8,[B(x)Bi(x) de. (12-9) 
j=0 


Comparison with the equations above shows that this approach is 
equivalent to selecting the weighting functions as 


wj(x) = B(x). (12-10) 


If the functions are complex, a better approach is to define the 
integrated squared error by the integral of the square of the 
magnitude of the error 


N N r* 
ISE = [ [ss -y £8(2)]| 50 -¥ 581) dx, (12-11) 
j=0 j=0 
from which the resulting equations can be shown to be 


N 
[BR x)e(x) de = Lg [BP (x)B(x)de. (12-12) 
j=0 


12.2. INDUCTIVE IRIS IN A WAVEGUIDE (REVISITED) 


We will now apply the method of mode matching to the same 
inductive iris in a waveguide problem considered before. The inci- 
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dent field is given by 
° _ WX 
E\(x, 2) = sin eta (12-13) 


where y, is 


" (12-14) 
The corresponding magnetic field is described by 
BiG 2) = gin (12-15 
rae Jon, a : ) 
The reflected (scattered) wave is described by 
= nrx 
E}(x,z)= Le, sin eve (12-16) 
n=1 a 
and 
= Ms nx 
Hy(x.z) =  —~e, sin em, (12-17) 
n=1 JORG a 


Because the total electric field is zero on the iris 
Ei(x,0) + E%(x,0)=0 for O<x<c, (12-18) 


which can be written 


= _ NWT Wx 
Ye, sin —=-sin— for 0<x<c, (1219) 
n=1 a a 


and because the tangential magnetic field is continuous where there 
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is no metal 
= Va _ ATX 
5 2, Sit —_—— = for ¢<x<a. (12-20) 
n=1J@Ko a 
Tf we terminate the summation at N, we get the approximations 
WV _ Nx _ Wx 
¥ e, sin —— = —sin— for O0<x<ec (12-21) 
P=) a a 
and 
Yn _ nwXx 
——¢, ai —— "0 for ¢<x<a. (12-22) 
n=1J@Rq a 


This is a special case of the problem discussed in Section 12.1, We 
could, for example, divide the interval (0, @) into N subintervals and 
use pulse functions for weighting, with N, pulse functions for the 
interval (0, c) and with N, = N — N, for the interval (¢, a). 


12.3. OTHER EXAMPLES OF WAVEGUIDE MODE MATCHING 


The concept of mode matching has been applied to a variety of 
waveguide problems. Any circuit that consists of a number of 
connected sections, cach of which is simple enough to have a 
solution involving known modes, usually can be analyzed by enfore- 
ing continuity of tangential field components at the interfaces 
between the sections. A change in width of a rectangular waveguide 
carrying a TE,, mode, as illustrated in Fig. 12-1, offers a simple 
example. The E and H fields on each side of the interface where 


ay 


ag 


FIGURE 12-1. Change in width of a waveguide. 
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the width changes easily can be expanded in infinite series involving 
well-known modal functions. Equating the tangential field compo- 
nents at the interface and terminating the summations to a finite 
number of terms yield several approximate equalities. Use of 
weighting functions, as discussed in Section 12.1, yields a set of 
linear algebraic equations. 


PROBLEMS 


1. For the waveguide width change illustrated in Fig. 12-1, write 
the appropriate modal expansions and indicate how a set of 
algebraic equations for the coefficients can be derived. 


2. For the situation in Fig. 12-1, assume an incident wave from the 
left. Using just the dominant modes and unity weighting. caleu- 
late the reflection coefficient. 


COMPUTER PROJECT 12-1 
For the situation shown in Fig, 12-1, with 


a4,=5cem 


a,=4cem 


assume an incident wave from the left which is a TE, mode at a 
frequency of 5 GHz. Write a computer program to determine the 
reflection coefficient. Use eight modes on each side, and match at 
the discontinuity using weighting functions which are sinusoids that 
correspond to the modes on the left side. Compare the result with 
the crude estimate obtained in Problem 2. 


Concluding Comments 


In this book we have described several basic methods of numerical 
analysis for solving electromagnetic problems. As can be seen from 
other books and the professional literature, many other methods 
have been derived. Two of these are the transmission line matrix 
(TLM) method—a time-domain method similar to the finite-dif- 
ference time-domain method and essentially a numerical version of 
Huygen’s principle—and the boundary-element method—a varia- 
tion of the finite-element method in which the boundary rather 
than the interior region is divided into clements. Additionally, many 
other new methods, most of which are combinations and variations 
on older methods, still are being derived. While much of this 
activity represents real progress, the variety of methods certainly is 
likely to be confusing. One question frequently aris: What is the 
best method for my particular problem?” The literature contains 
little in the way of real comparisons of methods, and what does 
appear is usually contradictory. The most important criteria in 
comparing methods are the time and effort to derive the algorithm 
and write the program, computer memory required, and computer 
run time required to achieve desired accuracy. This last criterion 
should receive more attention in the future. It seems to me that if 
there were really a superior method and everyone were rational, 
this method would be the only one used and the others would die 
away. Clearly this has not happened and hence a logical conclusion 
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would seem to be that different problems are best solved by 
different methods, 

Although computational electromagnetics has a fairly long his- 
tory, it still is a very active area of research and undoubtedly will 
continue to be. Most papers on numerical methods for electromag- 
netics are found in the three JEEE Transactions on Antennas and 
Propagation, Microwave Theory and Techniques, and Magnetics, and 
the interested reader is encouraged to read these journals and 
others, such as Radio Science. Improvements in computer hardware 
and software have influenced the development of the subject and 
expected improvements, especially those involving parallel compu- 
tation, promise an exciting future. 
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